Climate and Weather Risk in Natural Resource Models

This work, consisting of three manuscripts, addresses natural resource management under risk due to variation in climate and weather. In three distinct but theoretically related applications, I quantify the role of natural resources in stabilizing economic outcomes. Manuscript 1 examines policies for controlling the risk of outbreaks of cyanobacteria, algae blooms, stimulated by excess nutrient loading transported by uncertain and possibly changing rainfall patterns for a drinking water watershed in Rhode Island. The manuscript is based on a paper co-authored with Jim Opaluch. Manuscript 2 addresses welfare outcomes of optimal, as well as open access, extraction of a groundwater resource when demand is variable based on precipitation conditions. The groundwater stock is treated as depletable both in height and spatial extent. The model is applied to the Ogallala aquifer in Kansas. The manuscript is based on a paper co-authored with Todd Guilfoos. Manuscript 3 models the market for natural gas used by the electric sector in New England as a function of weather and pipeline capacity. The effect of changes in natural gas import capacity on natural gas price and quantity consumed as a function of the severity of winter weather conditions is quantified for the New England region. In Manuscript 1, we address policy designed to effect the risk of cyanobacteria blooms in a drinking water reservoir through watershed wide policy. Combining a hydrologic and economic model for a watershed in Rhode Island, we solve for the efficient allocation of best management practices (BMPs) on livestock pastures to meet a monthly risk-based as well as mean-based water quality objective. In order to solve for the efficient allocations of nutrient control effort, we optimize a probabilistically constrained integer-programming problem representing the choices made on each farm and the resultant conditions that support cyanobacteria blooms. In doing so, we employ a genetic algorithm (GA). We hypothesize that management based on controlling the upper tail of the probability distribution of phosphorus loading implies different efficient management actions as compared to controlling mean loading. We find a shift to more intense effort on fewer acres when a probabilistic objective is specified with cost savings of meeting risk levels of up to 25% over mean loading based policies. Additionally, we illustrate the relative cost effectiveness of various policies designed to meet this risk-based objective. Rainfall and the subsequent overland runoff is the source of transportation of nutrients to a receiving water body, with larger amounts of phosphorus moving in more intense rainfall events. We highlight the importance of this transportation mechanism by comparing policies under climate change scenarios, where the intensity of rainfall is projected to increase and the time series process of rainfall to change. The climate change scenarios show a shift towards a heightened risk of conditions supporting blooms and an increasing importance of spatial prioritization of nutrient control effort. In Manuscript 2, we introduce a new economic groundwater model that incorporates the gradual shift from irrigation to dryland farming as parts of an aquifer run dry. We accomplish this using an upside down cone to represent the spatial depletion, where the area of irrigable land above the aquifer shrinks as the water level decreases. Depletion of the aquifer may interact with uncertainty of the supply of water because the buffer that groundwater provides is no longer available. In this work, we identify the impact of spatial depletion on welfare gains from optimal management when rainfall is stochastic and follows a Markov process. Using a stylized model and dynamic programming, we estimate gains from moving away from current myopic extraction behavior to optimal use of the resource. When applied to Kansas over a section of the Ogallala Aquifer, we find gains from management ranging from 2.88% to 3.01% with larger gains achieved under uncertainty in the rainfall process. We find that including the dynamic of the gradual spatial depletion of the aquifer does materially impact welfare results compared to other estimates of the same region. Surprisingly the serial correlation of rainfall matters little. Empirically, multi-year droughts combined with the loss of access to the aquifer only slightly increases welfare gains due to the availability of dryland farming and the productivity of that option as a backstop when available. Manuscript 3 empirically estimates the effect of an increase in natural gas pipeline capacity in New England on monthly equilibrium natural gas prices and quantities for the electric sector. Weather plays an important role in defining the demand for natural gas due to its use for heating and electricity generation in the winter and through electricity demand for cooling in the summer. The cost of natural gas has important consequences to the wellbeing and cost of living for millions of customers either relying directly on natural gas for heating, or electric energy consumers indirectly. This paper presents results of reduced form price and quantity time series regressions using Generalized Least Squares (GLS) followed by results of a dynamic simultaneous equation model (SEM) of the market system. Using derived empirical relationships, prices and quantities for natural gas are estimated under various weather scenarios as well as under current and expanded capacity. I highlight the role capacity has in effecting the variability of the price of energy to the region. This work adds to the literature by providing empirical evidence and the quantification of the effect of constrained pipeline supply in an important energy market, where weather conditions, multiple demand sectors and alternative fuels determine the cost of energy. I find that capacity is a significant factor in the prices and quantities of natural gas consumed by the electric sector, with an increase in pipeline capacity of 1% leading to an average decrease in price of .48% and an increase in consumption of .2%. The SEM model finds both supply and demand to be price inelastic.

extent. The model is applied to the Ogallala aquifer in Kansas. The manuscript is based on a paper co-authored with Todd Guilfoos. Manuscript 3 models the market for natural gas used by the electric sector in New England as a function of weather and pipeline capacity. The effect of changes in natural gas import capacity on natural gas price and quantity consumed as a function of the severity of winter weather conditions is quantified for the New England region.
In Manuscript 1, we address policy designed to effect the risk of cyanobacteria blooms in a drinking water reservoir through watershed wide policy. Combining a hydrologic and economic model for a watershed in Rhode Island, we solve for the efficient allocation of best management practices (BMPs) on livestock pastures to meet a monthly risk-based as well as mean-based water quality objective. In order to solve for the efficient allocations of nutrient control effort, we optimize a probabilistically constrained integer-programming problem representing the choices made on each farm and the resultant conditions that support cyanobacteria blooms. In doing so, we employ a genetic algorithm (GA). We hypothesize that management based on controlling the upper tail of the probability distribution of phosphorus loading implies different efficient management actions as compared to controlling mean loading. We find a shift to more intense effort on fewer acres when a probabilistic objective is specified with cost savings of meeting risk levels of up to 25% over mean loading based policies. Additionally, we illustrate the relative cost effectiveness of various policies designed to meet this risk-based objective.
Rainfall and the subsequent overland runoff is the source of transportation of nutrients to a receiving water body, with larger amounts of phosphorus moving in more intense rainfall events. We highlight the importance of this transportation mechanism by comparing policies under climate change scenarios, where the intensity of rainfall is projected to increase and the time series process of rainfall to change. The climate change scenarios show a shift towards a heightened risk of conditions supporting blooms and an increasing importance of spatial prioritization of nutrient control effort.
In Manuscript 2, we introduce a new economic groundwater model that incorporates the gradual shift from irrigation to dryland farming as parts of an aquifer run dry. We accomplish this using an upside down cone to represent the spatial depletion, where the area of irrigable land above the aquifer shrinks as the water level decreases.
Depletion of the aquifer may interact with uncertainty of the supply of water because the buffer that groundwater provides is no longer available. In this work, we identify the impact of spatial depletion on welfare gains from optimal management when rainfall is stochastic and follows a Markov process. Using a stylized model and dynamic programming, we estimate gains from moving away from current myopic extraction behavior to optimal use of the resource.
When applied to Kansas over a section of the Ogallala Aquifer, we find gains from management ranging from 2.88% to 3.01% with larger gains achieved under uncertainty in the rainfall process. We find that including the dynamic of the gradual spatial depletion of the aquifer does materially impact welfare results compared to other estimates of the same region. Surprisingly the serial correlation of rainfall matters little. Empirically, multi-year droughts combined with the loss of access to the aquifer only slightly increases welfare gains due to the availability of dryland farming and the productivity of that option as a backstop when available.
Manuscript 3 empirically estimates the effect of an increase in natural gas pipeline capacity in New England on monthly equilibrium natural gas prices and quantities for the electric sector. Weather plays an important role in defining the demand for natural gas due to its use for heating and electricity generation in the winter and through electricity demand for cooling in the summer. The cost of natural gas has important consequences to the wellbeing and cost of living for millions of customers either relying directly on natural gas for heating, or electric energy consumers indirectly.
This paper presents results of reduced form price and quantity time series regressions using Generalized Least Squares (GLS) followed by results of a dynamic simultaneous equation model (SEM) of the market system. Using derived empirical relationships, prices and quantities for natural gas are estimated under various weather scenarios as well as under current and expanded capacity. I highlight the role capacity has in effecting the variability of the price of energy to the region.
This work adds to the literature by providing empirical evidence and the quantification of the effect of constrained pipeline supply in an important energy market, where weather conditions, multiple demand sectors and alternative fuels determine the cost of energy. I find that capacity is a significant factor in the prices and quantities of natural gas consumed by the electric sector, with an increase in pipeline capacity of 1% leading to an average decrease in price of .48% and an increase in consumption of .2%. The SEM model finds both supply and demand to be price inelastic.
vi ACKNOWLEDGMENTS I have many people to acknowledge that have taught, advised and supported me in a number of different ways. Todd Guilfoos has been a thoughtful, patient and insightful advisor over the past years. It has been a pleasure working with him on issues that we both find important and interesting, in addition to his guidance and perspective on making my way in the economics discipline. Through working with him as a research assistant, I am deeply grateful for the freedom he provided to explore problems as well as giving me the needed push to get ideas moving towards publishable work.
I would like to thank Emi Uchida for taking me in as a fresh graduate student and involving me in interesting and multi-disciplinary work early in my time at URI. I will forever be grateful for the opportunity to work with her and the rest international team in Tanzania. The experience has given me perspective and appreciation for how economics can be applied where natural resources and society are so intimately related. It was nothing short of life changing.

Introduction
Water quality impairments often present in discrete events, for example: fish kills, algal blooms, red tides, and beach closures. For the cases involving excess nutrient loading, these perceivable damaging occurrences are the more noticeable result of complex interactions across the landscape between humans, climate and ecological systems. Rainfall and subsequent runoff is the transportation mechanism for nutrients from the landscape to water bodies. Since rainfall is variable, any sought after water quality benefits of policy interventions are stochastic. Changes in policy create spatial and temporal configurations of human activity that in turn affect the distribution of conditions in a water body (Milon 1987). This paper assesses policies for controlling the risk of cyanobacteria (blue-green algae) blooms stimulated by excess nutrients. We examine cost-effective strategies that are based on addressing the probability mass in the upper tail of the distribution of phosphorus concentrations in a water body, which are disproportionately influential in supporting episodic blooms. Combining a hydrologic and economic model for a drinking water reservoir watershed in Rhode Island, we solve for the efficient allocation of best management practices (BMPs) on livestock pastures to meet a monthly risk-based as well as a mean-based water quality objective. We hypothesize that management based on controlling the upper tail of the probability distribution of phosphorus (P) concentrations in the reservoir implies different efficient management actions as compared to controlling mean conditions. We find a shift to more intense effort on fewer acres when a risk-based objective is specified with cost savings of meeting risk levels of up to 25% over mean-based policies. Additionally, we illustrate the relative cost effectiveness of various policies designed to meet this risk-based objective.
Rainfall and the subsequent overland runoff is a source of transportation of nutrients to a water body, with larger amounts of P moving in more intense rainfall events. We highlight the importance of this transportation mechanism by comparing policies under climate change scenarios, where the intensity of rainfall is projected to increase and the time series process of rainfall is projected to change for Rhode Island. The climate change scenarios show a shift towards a heightened risk of conditions supporting blooms and an increasing importance of spatial prioritization of nutrient control effort.

Background
Cyanobacteria blooms are challenging to address because they are episodic. They are patchy and difficult to predict (Hudnell 2008) because they depend upon the confluence of factors including sunlight, nutrients (especially phosphorus), pH, precipitation, water temperature, water flow, and water column stability (EPA 2012).
Under favorable conditions, cyanobacteria populations can rapidly proliferate and create extensive blooms dominated by a single (or a few) species that release toxins harmful to ecosystems and humans. Even in the absence of toxic effects, cyanobacteria blooms can result in adverse effects on taste and odor of drinking water (EPA 2012).
Of the conditions supporting blooms in fresh water, P loads and resultant in-lake concentrations are the most dependent upon anthropogenic sources and most amenable to management. P loads are stochastic since transport to receiving waters is mediated primarily by eroded soil moving in surface runoff (James & Alexander 1998; AL Heathwaite 2000; Morgan 2001). As a consequence, extreme precipitation events are influential in determining transport of P to water bodies, and contributing to conditions favorable to blooms. Policy within this stochastic context must be effective in controlling P loading in larger precipitation events that lie in the tails of the probability distribution. This could be of increasing concern in coming decades, as there are projections that climate change will increase high intensity precipitation events (IPCC 2007;Meehl et al. 2005;IPCC 2011;Kunkel et al. 2013).
The location and mix of nutrient control effort across a watershed affects both the average emissions as well as the shape of the distribution of emissions throughout varying weather conditions. When damages from P emissions are attributed to threshold based events, those damages are not a linear function of P concentrations, or continuous estimates of damages are unknown (Lathrop et al. 1998), then optimal nutrient control efforts that address mean conditions may differ greatly from those that address damaging occurrences or threshold violations (Shortle & Horan 2001).
Therefore, simple rules for improving abatement efficiency by focusing on the lowest marginal cost abaters in terms of average emissions do not necessarily hold (Shortle 1990). Rather, policy must consider the upper tail of the distribution of emissions and correlations across emitters.
Policy prescriptions to address a non-point source nutrient control problem in both a static as well as dynamic setting focus on location (and possibly time) differentiated taxes, subsidies, or standards based on estimated marginal damage of runoff, ambient conditions or input factors to approach optimal configurations of production and abatement effort (Goetz & Zilberman 2000;Segerson 1988;Shortle & Horan 2013).
Many studies use mean Nitrogen and Phosphorus concentrations or reduction targets as constraints to solve for the optimal allocation of nutrient control effort (Rabotyagov et al. 2010;Gitau et al. 2004;Arabi et al. 2008;Kling 2011). A number have explored costs in addressing probabilistic objectives of nutrient concentrations (Kampas & White 2004;Elofsson 2003;Shortle & Horan 2001;Gren et al. 2002;Milon 1987;Kataria et al. 2010;Zhu et al. 1994;Huang et al. 2012) and find distributional assumptions of emissions and the strictness of the probabilistic objective effect the efficiency of nutrient control design.
This work contributes to the literature by assessing the efficiency of the design of water quality improving policy to meet a risk-based objective. In contrast to previous work, we quantify the efficiency of policy designed to address a probabilistic constraint representing the risk-based objective as compared to mean-based policy using a hydrologic and numerical optimization model, while maintaining the important discrete nature of the problem. The resultant differences between efficient and second best allocations to address a risk-based objective show not just the total cost savings, but also the practical implications for how the allocation would differ across the watershed. By integrating climate model projections, we show that climate change may shift the risk profile and increase the importance of addressing the risk-based goal directly. The results underscore the need of spatially prioritizing nutrient control activities.

Conceptual Model
We assume a watershed partly consisting of heterogeneous agricultural fields with varying physical properties (e.g., slope, soil type, location). P accumulates in the soil over time and is transported off fields into waterways with precipitation events. P flows with water downstream to a water body, henceforth, reservoir, which provides ecosystem services valued by society.
We formulate the problem as finding the least-cost allocations of nutrient control effort to meet probabilistic constraints. This is in contrast to finding an economically optimal solution, which would include a damage function explicitly into the objective function. A damage function might take many forms depending on the water quality impairment. For instance, if damages were linear in concentrations, addressing mean concentrations would be sufficient. However, if the damages were convex to concentrations, or based on a threshold where damages are minimal and increase after a particular point, then higher order moments in the distribution of conditions matter to policy design. A probabilistic constraint can handle this type of least-cost problem by addressing the shape of the distribution of conditions (Shortle 1990). Our application to cyanobacteria blooms, which are known to be particularly reliant on high levels of P concentrations, fits this type problem. The damages are uncertain, but damaging blooms are, at the least, nonlinearly related to nutrient concentrations.
The distribution of total P loading at a water body is the sum of individual sources of emission less loss in transport. Stricter correlation between the sources leads to more burdensome probabilistic constraints and higher costs (Kampas & White 2004;Elofsson 2003). Given that each farm's emissions are highly correlated when rainfall is assumed consistent across the space, highly correlated loadings seem like a reasonable assumption, especially in a geographically small watershed. Each heterogeneous source will contribute a different but correlated distribution of nutrient loading based on its soil type, slope, inputs and location, or more generally a source's place on the P transfer continuum (Haygarth et al. 2005). We assume the watershed of concern is sufficiently small that rainfall is perfectly correlated across fields.
Assumptions of the distribution of pollutant load and transport can have large implications for estimates of efficient abatement (Gren et al. 2002;Kataria et al. 2010). In the numerical application below, we rely on a physical process-based hydrologic model to estimate the important relationships.
The goal of policy is to control the risk of cyanobacteria blooms by taking actions across the watershed to maintain P concentrations in a water body, a reservoir in our application, below a target concentration to avoid conditions that are favorable for blooms. However, P concentrations are stochastic, since loads depend upon random precipitation, especially upon extreme precipitation events. As a consequence, we specify policies that control the probability that a threshold is violated. Hence, the optimization problem is formulated as finding the least-costly set of control policies such that the probability of exceeding a threshold for P concentrations in the reservoir is less than a target level, .
Soil P concentration on a field at time , , is a state variable that evolves over time following the dynamic equation where F represents fertilizer inputs, A represents density of animal units, and represent the associated P application rates per unit, represents the loss rate for phosphorus (e.g., loss through plant uptake).
represents the field-specific rate of P emissions associated with a precipitation rate, r t, and soil P concentrations, , at time t. Therefore, represents the rate of P loss associated with soil erosion and runoff from rainfall. It is commonly assumed that the relationship between soil erosion, a major source of P loading, and rainfall intensity follows a power law. This assumption is based on the characteristics of sediment movement as a result of rainfall kinetic energy (van Dijk et al. 2002), and motivates the non-linear relationship between rainfall intensity and P loading. However, erosion and nutrient loss are a function of a number of other factors including soil saturation, slope and vegetative cover. In the numerical application below, we rely on a physical process-based hydrologic model to estimate the important relationships.
P control practices (BMPs) are represented as affecting runoff and are represented by the function . The simplest, and widely used formulation for P control is to use BMP efficiencies, which are a fixed percentage reduction of nutrient runoff based on edge of field studies. This is the assumption made in the empirical model presented below, but we will leave the theoretical model more general. Therefore, the P emissions, , in total loading (kg) terms from one field in time t are: where is the size of field i. Since the objective of the study is the water quality at the reservoir, the total loading is the sum of the loading from the n heterogeneous fields, as well as the additional loading from all other land use types that are not in the choice set of nutrient control. Because water remains in the reservoir for a period of time (residence time), the reservoir is itself a sink of nutrients and has its own in-lake processes. The P concentration, , in the reservoir at any one point in time is a stock variable progressing by an equation of motion.
is a function of both aggregate loading as well as the associated water volume, in order to convert the P loads (in kg) to concentrations (ppb).
represents the emission from all other land not considered for nutrient control actions.
The function represents the in-lake processes that might affect the progress of P concentrations.
Therefore, the stochastic optimization problem can be described as follows: Where: -Combination of BMPs on field i in time t -Set of BMP options -Cost of BMP employed on field i at time t -Size of field i -P soil concentration on field i in time t -Field specific rainfall to emissions function -Fertilizer applied on field i in time t -Animal density on field i in time t -P uptake and stabilization on field -P concentration in reservoir -P concentration threshold -P runoff of field i in time t -In-lake functions -P loading from land uses not in the choice space -Probability of exceeding threshold (risk level) -Rainfall represents the set of chosen controls given the available options (BMP), with costs in a given period . The choice of controls effects P runoff, , from each heterogeneous pasture land i with environmental factors (e.g., field size, location, slope, soil type) whose runoff to rainfall relationship is defined by . Each distribution varies by environmental factors and its shape is determined by choices of P control on each farm. The P threshold, , is the concentration of P that supports damaging cyanobacteria blooms. The term , represents the probabilistic constraint on the reservoir concentrations being above this threshold only % of the time or less. captures the risk-based water quality objective.
Replacing the probabilistic constraint above with the following represents the meanbased water quality objective: where captures strictness of the mean-based water quality objective.

Empirical Application
The framework above is applied to the problem of controlling conditions supporting A water body with a concentration exceeding 50 ppb, or 60 on the TSI scale, is at a level where algae blooms occur and where cyanobacteria may dominate. This would correspond to a yearly loading of roughly 1000 kg of P. However, this yearly concentration is made up of a series of monthly loadings and water flows. Therefore to construct the monthly risk measure for this paper and with the relatively short residence time of the reservoir, we use a monthly (May-October) inflow concentration above this 50 ppb threshold as the proxy for conditions that could support cyanobacteria blooms. This inflow concentration is calculated by using the variable monthly flow into the reservoir and the monthly total P loadings. The probability that the threshold is violated is the risk-level, in the conceptual model. This threshold and the cumulative distribution function of monthly concentrations modeled for the Regulating Reservoir can be found in Figure 2.
We use a flow-calibrated model, the Soil and Water Assessment Tool (SWAT), to estimate phosphorus loading to the reservoir as a function of rainfall and land use (Gassman et al. 2007). SWAT is a process-based water-modeling tool designed to simulate the surface and ground water quality and quantity, and to predict effects of land use change, land management practices, and climate change in complex watersheds. SWAT is a public domain model available free for download (http://swat.tamu.edu/). The SWAT model was developed and refined by the USDA Agricultural Research Service and Texas A&M University over a period of more than 30 years and has supported hundreds of scientific papers.
As part of a larger ecosystem service study in the watershed, our work focuses on practices installed on small acreage livestock farms to address nutrient loading and water quality in the watershed. Including variability in the animal density would be preferable as this information also describes the rainfall to P runoff relationships.
However due to privacy concerns, we were not able to obtain site-specific animal counts. As a consequence, we use the average animal density in the watershed (.5 AU/acre), as described by a livestock tally for the towns overlapping the watershed conducted by the Northern Rhode Island Conservation District (NRICD). This density was applied to all pastureland in the model. Given these assumptions, loadings from pasture land (about 3% of the watershed area) represent about 7% of the total P loadings to the Regulating reservoir.
In practice, nutrient control policies can only target a subset of loading to a given water body. There are natural loadings as well as other anthropogenic sources that may not be considered for management options. Therefore, policy in any watershed affects a subset of the total loadings and may have a limited effect on water quality.
However, certain loading sources may be more sensitive to rainfall intensity and runoff and nutrient control effort on this land may affect the distribution of conditions in a water body. Pasture land in particular is a large source of sediment bound P loading and a natural place to direct effort (Harmel et al. 2006 With three practices, this results in a total of 8 combinations on each farm. Each practice was characterized by a fixed and variable (per acreage) cost. The manure management BMP was translated from volume based variable cost to acreage based on average animal density in the watershed. These costs are shown in Table 1. For the simulation, we assume that it is a one-time decision to install and maintain the BMPs over their lifetimes.

Optimization
The objective of the numerical optimization is to solve for the least-cost allocation of integer optimization routine combined with output from SWAT and the cost and BMP effectiveness estimates described above to estimate the least-cost frontier.
GAs belong to a family of evolutionary algorithms, where possible solutions (individuals) to the objective function progress through populations and generations based on their relative fitness scores representing the optimization goal. Each new generation is created by crossing over (randomly combining) elements of better-fit individuals to create offspring. Through repetition, mutation and natural selection, the population converges to a solution. GAs have been used in nutrient abatement optimization frameworks by many authors (Gitau et al. 2004;Rabotyagov et al. 2010;Jha et al. 2009;Srivastava et al. 2002;Arabi et al. 2008). In our case, the solution is the allocation of BMPs across the landscape and its corresponding costs and risk level.
One of the inequality constraints is the risk level on concentrations of inflow to the reservoir. A GA is well suited for this type of spatial problem, as it does not require linearity, continuity, or differentiability in the objective or constraint functions (Arabi et al. 2008) The makeup of the initial population affects the efficiency of the routine as well as the width of the solution space that it searches, which is important in confirming a found solution is a global and not local optima. Taking advantage of some prior knowledge, the initial population for the GA was seeded with individuals representing uniform adoption of each practice or across the board adoption of all possible BMP combinations, as well as allocations found by incentive based policies described below. The remaining individuals were generated randomly using a uniform distribution. In order to cover a wide range of solutions, the algorithm was run with a population of 500 individuals, for a maximum of 500 generations, or until the variability of the fitness values of individuals within the population was within a function tolerance.
We solved for the least-cost allocation of BMPs for each level of risk within the lower and upper bounds of our choice set, which are the risk levels between no installations and all practices installed on every farm. We proceed by seeding the initial populations for the GA with previously found optimal solutions. We execute this process under the probabilistic constrain as well as under the mean constraint for each rainfall scenario. As a robustness check and to confirm the GA was behaving correctly, we compared the allocations found by the GA under the mean concentration constraint (without the allocations resulting from policy scenarios in the initial population) to the simple method of sorting the farms by the average per unit cost of abatement, which defines the efficient progression for the mean concentration goal.
The two solutions were identical.

Policy Scenarios
We construct cost functions for each of the following policy designs, where the function for each design indicates the cost of achieving each risk level. We compare the efficient allocation found using the optimization process above to three simulated policies: a location differentiated subsidy; a flat subsidy based on the relative effectiveness of the BMPs, but devoid of location specific factors; and a uniform standards approach. The subsidies are simulated assuming each farmer chooses the practice(s) on his/her farm to maximize the value of a one-time up front subsidy minus the discounted present cost of installing and operating the practices. We assess efficiency in terms of the cost of the installed practices to meet risk levels. The subsidy is treated as a transfer payment rather than a social cost.
The location-specific subsidy assumes a one time up-front payment is offered for installing practices, with the total payment based on the average annual amount of P abated at the reservoir due to the practice or combination of practices on a farm. The payment is fixed across the watershed for a kg of P abated, but the abatement effects of BMPs vary for a particular practice on a farm. For example, if farmer x installs a practice which would result in a average reduction of 1 kg P/year from reaching the reservoir and the watershed wide per kg P incentive is $1000, the farmer weighs the total one time payment of $1000 against the cost of installing and operating that practice over 20 years. 2 The farmer chooses the combination of practices that maximizes net profit, including the option of no BMP installations. Clearly, this type of subsidy requires an agency to model and offer site-specific subsidy schedules based on a particular farm's heterogeneous effect, which is a tall task, both in terms of transparency, politics and operational constraints. The benefits, quantified in our results, would be a more cost-effective allocation of effort due to using the sitespecific factors.
This setup is similar to the USDA's Conservation Reserve Program (CRP), where land is retired and put into environmentally friendly management practices in return for a rental rate from the government. The program relies on bids from the farmers that are then ranked by an Environmental Benefits Index (EBI). The EBI balances many environmental goals such as soil conservation, nutrient retention and wildlife habitat.
The EBI is spatially differentiated by distance to water bodies, soil type and slope, as well as giving weight to certain land based on particular initiatives. The government picks the top ranked bids based on cost-effectiveness and budgetary constraints. If accepted, the yearly payments are set at the bid price the farmer submitted for the life of the contract (USDA 2015). The difference compared to the first location-specific subsidy policy is that we assume a payment is offered based on the farm's effect and the farmers choose if that is in their best interest. The equivalent in the CRP framework would be if the government offered a fixed payment based on a unit of EBI and informed a farmer of their possible payment for installing a practice on their 2 An alternative way to think about this on time upfront payment would be to consider an equivalent yearly payments whose present value is equal to this $1000. particular farm based on that farm's total EBI. Thus, there is a difference in market mechanism, with the CRP program based on accepting asking prices from farmers, and the subsidy described in this paper based on farmers accepting bid prices from the government. These differences could have implications to participation, size of transfer payments, but not necessarily allocate efficiency in terms of social cost.
The second policy scenario, a flat subsidy, assumes a subsidy based only on information of farm size and practice, but no other site-specific information.
Therefore, it is a subsidy based on an assumed per acre relative effectiveness of the BMPs and the farm size. Hence, if the payment for a practice were $100 per acre, a practice that is twice as effective in general would receive $200 per acre regardless of the particular farm's P loading. The installation decisions follow the same logic as above, with the farmers choosing the practices that maximize their net profit.
These subsidy based policy designs are simulated for the full range of achievable allocations by varying the payments to cover the whole range of effects. This would be the per kg P payments that would induce the range from the minimum to where maximum BMPs are installed based on the implementation choices by farmers. The allocations are evaluated based on achieving risk levels at various costs in the results section.
The final sets of policies are uniform standards. Under the uniform standard policy, each farm is required to adopt a particular combination of the three BMPs. A cost function for uniform standards is constructed by calculating the costs and associated risk levels for each of the eight possible BMP combinations.

Climate Scenarios
As indicated above, extremes in precipitation intensity are central for determining P loads. This could have important implications for P concentrations and associated risk of cyanobacteria blooms because, while significant uncertainty is involved, there is an expectation that the Northeastern US will experience an increase in total annual rainfall and an increased percent that falls in extreme precipitation events (Karl et al. 2009;IPCC 2011). We use the process described below to model rainfall parameters We model precipitation with a two-stage process. First, we determine whether a day is wet or dry. Then, we determine precipitation amounts for wet days. The probability of a wet day (more than trace amounts of precipitation) is modeled using a first-order Markov process (Srikanthan & McMahon 1999). Hence, the first stage is to model the probabilities of rainfall on date t for the cases where date t-1 was wet and where date t-1 is dry. To do so, we specify the conditional probability of date t being wet when date t-1 was wet (P(w t |w t-1 )) and conditional probability of date t being wet when date t-1 was dry (P(w t |d t-1 )). We model rainfall amounts for wet days using a two-parameter gamma distribution following the methods described in Davison et al, (2005) and Geng et al. (1986). We estimate the parameters of the gamma distribution of rainfall applying method of moments matching to historic daily rainfall data and the climate model runs.
We use historic daily rainfall data for Rhode Island from 1961 to 2010 from the Kingston, RI weather station to calculate these two conditional probabilities and gamma parameters for each month for the baseline rainfall scenario. We then create a 100-year time series of simulated daily rainfall with the observed statistical properties for input to the hydrologic model. We use a longer series for the simulation than the 50-year windows used to calculate rainfall parameters in order to better populate the discrete distribution of monthly concentrations under each scenario. To create the future rainfall scenarios, we calculated the parameters described above for both historic   These results are summarized in Table 2. The average of the changes in the models imply that the intensity of daily rainfall will increase, which can be seen in the percent increase in the scale parameters from the climate change models. This shifts the rainy day distribution to increase the frequency of larger rainfall amounts. The decreases in the shape parameter implies rainfall on a wet day will be less normally distributed and become closer to an exponentially distributed random variable, in other words: less symmetric and more skewed with more probability mass in the upper tails. These parameter changes mean an increase in the intensity of daily rainfall events. The change in the Markov chain parameters are mixed, but a decrease in the probability of 4 We acknowledge the World Climate Research Programme's Working Group on Coupled Modeling, which is responsible for CMIP, and we thank the climate modeling groups for producing and making available their model output. For CMIP the US Department of Energy's Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. 5 Radiative forcing is the amount of energy that is retained by the atmosphere as a result of greenhouse gases. A great explanation is found in Chandler (2010) a wet day followed by a wet day, P (w t |w t-1 ), means that when rainfall occurs it will occur in shorter more intense events over a shorter run of consecutive days than under current rainfall conditions. The probability of a wet day following a dry day, P (w t |d t-1 ), increases, which means dry spells will be shorter with the exception of August where the probability slightly decreases. In summary, these climate models point to a climate that is wetter in total and with more intense rainfall events than under the current conditions. These more intense events should cause more overland flow, sediment erosion, leading to nutrient loading at the reservoir. We rely on the SWAT model to simulate the time series nature of soil saturation, runoff at rainfall events, as well as nutrients transported and associated water flow to the reservoir.

Results
The simulation and optimization provide the least-cost allocation of BMPs and those induced by the policy options in the watershed as well as the resultant distribution of monthly total P concentrations at the reservoir. We compare the allocations in terms of cost, density of practices and size of participating farms. The following results are broken down into those under the observed rainfall scenario followed by those under the climate change rainfall scenarios.

Baseline Rainfall Scenario
The range of outcomes as a result of BMPs installed on pasture land in the watershed is in Figure 2. The minimum BMP adoption line (thin black line) shows the reservoir's P concentration exceeding the 50 ppb 18% of the summer months (May-October), while at full BMP adoption (thick black line) exceeding this threshold 12% of summer months. The distance between the two lines in Figure 2 at the 50 ppb threshold represents the whole range of feasible outcomes as a result of pasture land BMP implementation in this watershed. The relative cost efficiency of the various policies differs across risk levels, which can be seen by the distance between the various alternative policies and the efficient frontier changing over the range. The allocations resulting from location-specific subsidy, taking advantage of locations' differentiated effects, are the closest to the efficient frontier inducing allocations similar to the least-cost allocations. The remaining difference represents the efficiency of addressing the risk-based objective explicitly.
Take an example of achieving a risk level of 13%. The cost of moving from the current state of 18% to 13% could be achieved optimally at $300,000, but at $335,000 if allocations were designed to meet mean concentration goals instead. This is a savings of roughly 12% by adopting allocations based on the risk-based objective. The potential cost savings of explicitly adopting a risk-based allocation ranges from 0 to 25%, with an average of 6.5%.
We have established that the efficient allocation of BMPs across the watershed differs in cost for the risk-based goal versus a mean-based goal. Interestingly, BMP adoption strategies for the two approaches differ. The efficient set of BMPs under the riskbased approach is more targeted, with more intensive controls placed on fewer acres.
By comparing the differences in allocations, we can show which farms and practices change at each risk level. Figure 4 shows the differences in total acres in each practice.
Noticeably, the optimal allocations have fewer acres in practices in total. In terms of total acres of installed practices, the risk-based objective results in more targeted effort, with fewer acres enrolled. This varies, but the number of acres averages 3% less while increasing to 15% less in the higher effort and lower risk levels. The biggest differences are in the nutrient management practice, a relatively cheap and less effective practice, with fewer enrolled acres under efficient allocations to meet a riskbased goal.
An additional way to see the differences is using a heat map approach (Wilkinson & Friendly 2009), where we characterize the discrete differences based on farm size between the second best and least-cost allocations. Figure 5 shows the discrete changes between the two allocations, sorted and shaded by their relative size. Each rectangle represents a binary decision for a BMP on a farm. The farms are sorted left to right by size from smallest to largest. The darker color implies a larger farm with a practice installed in the risk-based allocation but not under the mean-based allocation.
The lighter color implies a larger farm without that BMP in the risk-based allocation, but existing in the mean-based allocation. The neutral color means the BMPs are the same in both allocations. There is some variation across the risk levels, but in general we see the shift to fewer farms in more expensive and impactful practices and a shift from nutrient management (less impactful, less expensive practice) to manure management (more impactful, more expensive practice). We also see many smaller farms shifted out of the risk-based allocations to BMPs on fewer larger farms. The lower the risk level, the more farms switch their practices, which means the differences are largest in allocations near the upper limit of what is achievable for risk reductions.
Given that P loading from pasture lands only made up 7% of the total modeled loading to the reservoir, the range of reduction in mean P concentrations from installing practices (23.9 to 27.3 ppb) seems small and the costs hard to justify depending on what damages might be associated with changes in mean concentrations. However, depending on the damage function, the effect of nutrient control efforts may look more promising in having an overall impact on the water quality by lowering the chances a water body supports damaging events. For example, over a 20-year timeframe, the range of effects in terms of risk levels means moving the reservoir from 21.6 summer months to 14.4 months with P concentrations above the 50 ppb threshold.
While we treated risk based on strict threshold violations, it may be the case that the shape of the distribution of conditions below the threshold matters for other water quality goals. It is informative to see the tradeoff of targeting the risk-based objective in terms of other moments of the distribution of conditions. There exists a tradeoff between addressing the risk levels efficiently versus the mean concentrations. The shift towards the addressing the risk levels most cost-effectively comes at the expense of a slightly higher mean concentration at each risk level when the allocations differ.
This difference is small and between 0, where the allocations match, and .03 ppb at the maximum difference in spatial allocations of BMPs. So, while there is a difference, shifting policy to address a risk-based goal is nearly costless in terms of increased mean concentrations.
We showed that the type of objective, mean or risk-based, matters to the cost-effective adoption of BMPs across farms and the magnitude of the efficiency of policy designs.
The risk-based objective leads to more intensive control on fewer acres, and a shift from more spread out effort to fewer, but larger farms with installed BMPs. Our project was focused on BMPs, where the cost structure had fixed and per acre variable cost components. As a result, the size of the farm to which the policy is aimed matters to the overall cost, with larger farms having generally a smaller per unit cost of P reduction all else equal. This would also be the case where significant economies of scales occur for pollution control measures. Economy of scales applies to the cost of controlling not simply the mean but also higher order moments of a farm's runoff distribution. This lead to a trend of targeted and concentrated effort on fewer but larger farms. Therefore, policy that is spread out more evenly, standards, for example, fared worse.  Comparing optimal allocations between baseline and the RCP 6.0 climate scenario, the allocations are similar in terms of the progression of practices installed and the relative amount of each through the risk levels. The climate scenarios do not imply a major difference in efficient policy or allocations, but imply larger gains to targeting the riskbased objective explicitly, with larger average gains and a larger maximum gain, 46%

Climate Scenario Results
versus 25%. As rainfall becomes more intense, the possible benefits from moving towards efficient policy based on a probabilistic objective become larger and uniform standards become progressively less desirable. Given the differences in allocations, with a shift to fewer but larger farms with BMPs installed when addressing a riskbased objective, the results of the climate change scenario suggests the gains from doing so are larger. Under future climate conditions, spatial prioritization may become more beneficial for this threshold based water quality objective.
We presented a middle of the road emissions pathway, RCP 6.0, above. From Table 2 we can see the changes under a more severe scenario of climate change. The same trends in the cost effectiveness are intensified under the RCP 8.5 scenario. Figure 8 shows the cost at risk levels for the baseline and two climate change scenarios together. The risk levels for the two climate scenarios overlap, so a direct comparison at risk levels can be made. For a given risk level, the allocations under the RCP 8.5 scenario are further up the cost curves and result in more BMP implementation.
Equivalent risk levels are achieved at costs $100-500 thousand more than under RCP 6.0 scenario, visually, the distance between the two curves in Figure 8.

Conclusion
The variability of natural systems poses challenges in addressing and quantifying the costs and benefits of policies. This paper uses a combined hydrologic and economic model to solve the discrete optimization problem to compare the cost effectiveness of policies for nutrient loading to the efficient frontier under both a risk-based and a mean-based water quality objective. Controlling a risk factor for cyanobacteria blooms in drinking water, which are dependent on conditions violating a eutrophic threshold, motivated the approach.
Achieving a risk-based water quality goal, especially applicable to cases where damages are convex to a pollutant or are based on a threshold, results in a different management strategy than meeting objectives based on average conditions. Different sets of pollution control actions and allocations across space can affect higher order moments of the distribution of conditions as opposed to just the mean. We find that achieving the risk-based objective leads to a more targeted approach, with more intensive measures being adopted on fewer but larger sites. Given the results above, focusing on fewer, but more important larger farms may be more cost effective in controlling episodic loadings and damages. This is in contrast to traditional approaches of uniform technology standards or spatially uniform incentives.
Policy designs vary in the amount of differentiation in terms of how they treat heterogeneous sources of emissions. For a risk-based objective, the importance of spatial prioritization becomes larger. In our application, the efficiency of three alternative incentive policies showed the relative cost effectiveness of various second best policies to the risk-based efficient allocations. Utilizing location-specific mean loadings, and a watershed wide uniform payment for P abatement, matched the optimal allocations to address a mean concentration goal. This was closest to the riskbased efficient frontier out of the alternatives considered but still up to 25% more costly to meet a given risk level than the efficient allocations.
Rainfall, which is the transportation mechanism for many pollutants affecting water quality, is projected to change in its intensity and timing with climate change. This will affect both the mean and higher order moments of the distribution of water born pollution conditions in the environment. We find this has implications to the efficiency of policy. Specifically, the climate change scenarios increase the cost difference between mean and risk-based policies. Incorporating location specific factors into policy when the damages are based on a threshold is more important under more intense future rainfall. For our application, this difference grows with up to 46% lost between the efficient and second best policies. The non-location differentiated policies and uniform standards are also progressively worse in terms of efficiency.
Depending on how often and how far away from the threshold a water body of interest is, marginal changes in effort can have effects on the risk profile of a water body. To best make use of what might be a limited area of land where nutrient control effort might take place, prioritizing locations based on addressing the objective that best represents the damage or goal of the policy could save scarce resources. For the case of cyanobacteria blooms, we find that policy designed to address the probability of crossing a monthly eutrophic threshold of phosphorus concentrations differ than that to address mean concentrations, resulting in more concentrated actions being taken on a smaller number of sites.

Acknowledgments
The authors gratefully acknowledge the funding for this project from U. S.   NOTE-Gamma parameters, shape and scale, as well as the Markov chain parameters of wet days were fit to observations from for 1950-1999 by month in mm of rainfall from the Kingston RI weather station. The differences between past  and future (2050-2099) runs of climate models (both RCP 6.0 and RCP 8.5 emission paths) for these same parameters are presented as adjustments. These % changes were applied to the observed parameters that were then used to create the time series of future rainfall for climate scenarios.

FIGURE 2 CDF of Monthly P Concentrations at the Reservoir
Baseline Rainfall Scenario NOTE-The x-axis represents monthly P concentrations of inflows to the reservoir. The y-axis is the probability that the reservoir is below that concentration. The two lines represent the outcome of the minimum and maximum BMP implementation on pasture land in the watershed. The vertical line is the 50 ppb threshold used as an indicator of a high risk of algae blooms. The distance between the two lines at the threshold is the feasible risk levels, which is the x-axis in Figure 3.

FIGURE 3 Total Costs at Risk Levels Baseline Rainfall Scenario
NOTE-The solid lines represent the cost of achieving risk levels assuming optimal policy (solid black), a payment based on the yearly average P reductions as a result of installed practices (dashed), a flat per acre payment based on the relative BMP efficiencies (solid grey), but not taking into account other site specific characteristics. The black dots represent the risk levels and costs assuming standards.

FIGURE 4 Differences in Acres with BMPs
Risk-Based Compared to Mean-Based NOTE-The differences in the total amount of acres in the three practices (prescribed grazing, manure management, nutrient management) are shown between the monthly risk and meanbased objective. The positive acres represent additional acres in practices for optimally addressing the risk-based objective. The zero line shows the two allocations of BMPs to be the same and negative implies more acres in practices than under the mean loading objective.

FIGURE 5 Discrete Differences in Installed BMPs
Risk Versus Mean Allocations NOTE-Discrete changes between the two allocations at various risk levels, sorted and shaded by their relative size. Each block represents a discrete decision on a farm. Darker color implies a larger farm added to the risk-based allocation and whiter color shows a larger farm subtracted from the risk-based objective, or alternatively thought of as existing in that practice under the mean-based allocation, but not under the risk-based allocations. There is a shift towards more concentrated BMPs on larger farms under the monthly risk-based objective. .13 .12 .14 .15 .16 .17 Risk Level Acres

Introduction
Groundwater plays an important role in mitigating the effect of weather variability on economic activity. Agriculture, in particular, is highly dependent on rainfall and is thus uniquely exposed to weather-related risk. Many groundwater aquifers around the world are stable supplies of water for irrigation to make up for deficits in natural rainfall. As more than an additional source of water, they act as a water buffer in low rainfall years. According to the World Bank, it is estimated that 2 billion people worldwide depend on groundwater for drinking water. In addition, 50% of cereal production depends on groundwater for irrigation (Wijnen et al. 2012 Irrigated agricultural plays an important role in reducing the lifespan of many of these aquifers. For instance, over the Ogallala aquifer, irrigation began in the late 1800's but intensified greatly after WWII with the introduction of pivot irrigation systems. As of 2011, total storage across the aquifer has fallen 8.3% since 1950. This depletion of the resource has not been spatially uniform across the whole system with some limited areas gaining while most areas have faced falling levels. In some places, decreases as much as 242 feet have been observed, as can be seen in Figure 1 (McGuire 2012). In western Kansas, areas are left with little or no saturated thickness, and farming is forced to transition from irrigated crops, mostly corn, to non-irrigated crops such as sorghum, wheat and cotton (Steward et al. 2013;The Economist 2013). This switch from high value, stable, irrigated crops to lower value dryland varieties is captured by the spatial stock dynamics of our model.
As groundwater is depleted, there is a loss in its marginal value due to higher costs of extraction. This is in addition to a loss of the groundwater stock's marginal impact in reducing the variability of returns to farming (Tsur & Tomasi 1991;Knapp & Olson 1995). Despite this buffer value being a potentially large percentage of the total value of a groundwater resource 6 , welfare gains estimated when moving from open access to optimal behavior are relatively small in size, supporting the findings from Gisser and Sanchez (1980). 7 A review of research addressing the curious finding of Gisser & Sanchez (1980) by Koundouri (2004)  6 Buffer values were estimated at between 5-84% of total value, but a comparison between myopic and optimal management was not presented (Tsur and Tomasi 1991) 7 Knapp and Olson (1995) find gains to optimal management of 2.6% for Kern County, CA programming model. They find that the scarcity value of water increases tenfold under drought conditions even when optimally managed, although they do not present results in net present value form. The method implies that the exact time series of rainfall is perfectly known over the optimal planning horizon, which a stronger assumption then knowledge of the rainfall process, but not the exact realization. Roseta-Palma & Xepapadeas (2004) present a more general formulation and use robust control techniques to estimate optimal groundwater resource use. Precautionary behavior implies that optimal management under uncertainty results in extraction rates that are lower to save for future shortcomings than under more certainty in the random process of rainfall. Welfare effects of optimal versus open access management are not discussed in this framework, and the results are theoretical and not parameterized to a groundwater source. Zeitouni (2004) presents results for an aquifer with stochastic recharge but constant, not rain dependent, demand and constant marginal benefit. This leads to a bang-bang solution for optimal control and a target groundwater level. In contrast, our paper will use constant recharge, crop yield functions and duel sources of water, rainfall and groundwater, in describing the marginal benefit curves, demand, extraction rules, and resulting stock dynamics.
In addition to uncertainty in rainfall and surface water sources, other papers have addressed irreversibility and catastrophic loss of a resource. Tsur & Zemel (2004) and Leizarowitz & Tsur (2012) model the threat of a discrete permanent change in the system both with certain and uncertain stock dynamics. In specifications of stockdependent event risk, prudence or lower extractions are optimal to save for future states, although this does not hold for irreversible exogenous events, which imply larger extractions than otherwise. In contrast, our paper does not address discrete event risk or regime changes. We will compare optimal policy and resource savings under specifications in the rainfall process as well as quantifying the implications to the externalities associated with open access usage. over, representation of the pumping height externality when assuming a uniform single cell aquifer, known as the bathtub model. The specification of this particular externality is meaningful for welfare analysis since it is the magnitude of the pumping height externality, when taken into account in optimal management, which results in welfare gains of optimal extraction. Depending on the spatial and size distribution of users, welfare gains to management can be over or under estimated when using a bathtub type model of groundwater. In addition, the degree of the associated externalities dissipated with increasing degrees of strategic behavior and ownership by extractors, as pointed out by Brozović et al. (2010) and Guilfoos et al. (2013). A different, and possibly larger stock externality, that of the depletion of the total spatial extent of an aquifer, is explored in this paper with a more detailed multi-user spatial approach put aside for clarity. Simply put, in the model presented, as a groundwater resource is depleted, so too is the irrigable land lying above with access to the aquifer.
Thus, farm operations on the land that have lost access to groundwater must transition from irrigated crops to dryland crops. 8 We introduce a new economic groundwater model that incorporates the gradual shift from irrigation to dryland farming as parts of the aquifer run dry. We accomplish this using an upside down cone to represent the spatial depletion, where the area of irrigable land above the aquifer shrinks as the water level decreases. Depletion of the aquifer may interact with uncertainty of the supply of water because the buffer that groundwater provided is no longer available. In this work, we identify the impact of spatial depletion on welfare gains from optimal management when rainfall is stochastic and of a Markov process. With this model, we estimate gains from moving away from current myopic extraction behavior to optimal use of the resource. When applied to Kansas over a section of the Ogallala Aquifer, we find gains from management ranging from 2.88% to 3.01% with larger gains achieved under uncertainty in the rainfall process. These results quantify the importance of the spatial externality as well as including uncertainty. In the process, we make two contributions to the literature. First, the introduction of a model that allows for a gradual spatial When dynamics of the loss of irrigable land are included, the buffer value of groundwater in reducing the variability of profits may be particularly important, as more farmland is transitioned to dryland, rain dependent farming. The extent to which groundwater can act as a buffer against variable future rainfall is also a function of the costs of extraction in a given period, which in this model, is dependent on the stock, or height of the groundwater. Therefore, with higher levels of groundwater, more can be extracted across a larger area of farmland to meet demand in realizations of low rainfall, thus limiting the variability of aquifer-wide profits in any given year. As groundwater levels fall, so too does the amount available for withdrawal in response to low rainfall. Therefore, groundwater depletion increases the variance of possible realizations of aquifer wide profits in subsequent periods. The extent of the welfare gain due to optimal management is an empirical question, which we address with an application.

Spatial Externality
We model the dynamics of how groundwater evolves and access to the resource changes by introducing a spatial depletion function where the amount of irrigable farmland above the aquifer is a function of the pumping height ( Figure 2). This specification represents an upside-down cone shape where the area above the remaining groundwater faces uniform pumping heights, but the area no longer above the groundwater is left with without any option for irrigation. Doing so creates an analytical relationship between the amount of water extracted from the aquifer and the evolution of the height of groundwater and the changing access to groundwater. This Where is the initial groundwater height, the radius of the circle that is the base of a cone that represents the irrigable land area above the aquifer. Derivation of equation (1) can be found in the Appendix. This specification captures both the heterogeneity in saturated thickness as well as allowing for changes at the extensive margin (switch to dryland practices). This model is most applicable to areas where the marginal cost of increased pumping heights is relatively small compared the value of groundwater, or when the spatial distribution of wells and water demand is relatively uniform. By simplifying this spatial externality, it allows us to investigate the stochastic dynamics while keeping the size of the stochastic dynamic programming problem computationally manageable.

Rainfall Process
There are many ways to model the variability of rainfall, depending on the time and spatial scale appropriate to the problem (Srikanthan & McMahon 1999). We choose to model annual rainfall expectations and realizations two ways: the first scenario assumes independent and identically distributed (i.i.d.) realizations of rainfall. There are discrete rainfall states, , and corresponding probabilities, . The probability of next year's rainfall, , is independent of the current rainfall state, , thus . (2) Equation (2) defines expectations when rainfall is i.i.d., where is the expectation operator.
However, it is clear that rainfall time series are not serially independent and there exists many stochastic rainfall generation methods to model the serial dependence (Srikanthan & McMahon 1999). Therefore, the second scenario presented in this paper assumes a simple Markov chain process in order to replicate climate persistence, particularly droughts, where the probabilities of future rainfall states are a function of the current year's rainfall. This process defines probabilities of the future state that are conditional on the current state, , and By specifying the transition probabilities of moving from one state to another, we capture the time dependent process and the persistence of annual rainfall.

Economic Benefits
In order to estimate the per-acre return to farming, we assumed a section that is irrigated acreage and a section of farmland that is dryland acreage. The returns to irrigated land are a function of per-acre yield, price, irrigation pumping cost and the quantity of irrigation water extracted and applied. Equation (4) gives the per-acre returns to irrigated farmland.
Where is the per-acre yield, is the irrigation water applied per-acre, is the price of the irrigated crop and is the lift, , dependent marginal cost of one acre foot of groundwater extraction. is the elevation of the groundwater level as compared the surface elevation, Taken together, the cost function gives the marginal cost of pumping one acre-foot of water to the surface for a given pumping height. Myopic farmers that irrigate do so to maximize profits in each time period by choosing to maximize (4). Farmland that does not have access to the aquifer only has rainfall as a water input into the production of crops and is represented by, Where is the per-acre yield of dryland crops as a function of rainfall and the prices of those crops.
As mentioned before, irrigable farmland is itself a function of pumping height. So a single period's aquifer wide return would be the area-weighted sum of irrigated and dryland profits: We assume the homogenous per-acre irrigation water demand across the irrigated portion of the aquifer, so the total water extracted is: Based on the cumulative extraction , the groundwater height changes through time following the equation of motion: R is the natural recharge, the percent of irrigation water returning to the aquifer, A the total area of the aquifer, S the storitivity and the percent of irrigated land over the initial aquifer area from equation (1) The social planner's problem is to find the extraction path to maximize the discounted sum of future profits, subject to the equation of motion (8) and physical constrains, as given here: s.t.

Where
, is the expected value operator over the stochastic rainfall variable, . is the discount factor, the year index and and are the maximum and minimum groundwater levels. The random nature of rainfall is important in this setting as groundwater is used to augment natural precipitation in dry years. either follows the deterministic, i.i.d scenario, or the Markov chain assumption highlighted above.

Dynamic Programming Problem
Since the objective of this paper is to implement stochastic rainfall to address welfare gains from optimal management, we use a dynamic programming approach (Bellman 1957). A discussion of dynamic programming as applied to a groundwater setting can be found in Provencher & Burt (1994). The value function represents the expected present value of future benefits of the system assuming optimal management in all future periods.
By applying the principle of dynamic programming, the first order conditions for this problem are given by the Hamilton-Jacobi-Bellman equation for this discrete time stochastic model (Brito 2007).
Where is the discount factor equal to , and the discount rate. is the expectation operator over the random variable and the value of next period's stock assuming optimal behavior in all subsequent periods.
represents the transition equation as a function of current extraction decisions and groundwater height, which in our case is the same as equation (8) above.
Intuitively, the optimal extraction of groundwater balances the marginal benefit today against the discounted marginal benefits in all subsequent time periods given the expectations of the random variable. The first order condition with respect to withdrawals for the Bellman equation above implies the following Inputting equations (6 and 8) into equation (12), Although we do not know explicitly, along the optimal path, the marginal value of extraction today (the first term) should equal the discounted marginal cost (the second term) imposed by period 's extraction. The envelope theorem implies that the derivative of the value function with respect to the groundwater stock level generally is: Solving for in (13) and substituting that into equation (14), moving ahead one period and rearranging terms, we arrive at the discrete time Euler equation, which does not contain the unknown value function.
This represents the tradeoff between withdrawals today, the first term, and the opportunity cost, the part after the equal sign. With our functions, The opportunity cost includes many components.  (Gisser andSanchez 1980, Feinermann andKnapp 1983). Using these restrictive assumptions gives us the equivalent to equation (16) for the bathtub model.
We can see the term , appears in equation (16) but not (17) and is positive when , implying an additional opportunity cost of extraction over the bathtub model. The term in equation (16), , represents the additional impact on future benefits from the change in the water table in the spatial depletion model. If this is greater than unity it provides an additional cost in the spatial depletion model compared to the bathtub model. It is unclear if this term is always greater than one and provides an addition cost to extraction compared to the bathtub model. It is clear that when is larger that there are greater costs than when is smaller given assumption about the rate of depletion. 9 However, a result of the cone assumption is that as the water height is lower in the aquifer, a similar recharge will result in a larger height change due to the reduced size of the volume it is filling. Therefore, the total effect of adding the spatial externality is a matter of parameterization. 9 The key assumption here would be that the second derivative of the depletion function is constant, and the rate of change in irrigable area were constant, then the costs of the change are higher when γ(x)is relatively large because it only enters in the numerator of equation (16).
Solving equation (11) means finding the function either explicitly or numerically. The addition of a stochastic element and a non-linear spatial depletion function makes analytic solution intractable in this case (Brito 2007 The physical parameters of the model are presented in Table 1. The spatial externality is modeled using the function in equation (1), which captures the loss of irrigated acreage as groundwater levels fall. The initial aquifer area is set at 2.19 Million acres, of which 373,200 acres are irrigated or 17% of the total land. We assume that this total acreage of farmland (irrigated plus dryland) remains constant as the aquifer area is depleted. A discussion of the implications of allowing the amount of farmland to vary is found in the discussion.
The bottom of the aquifer is set at 2,892 feet above sea level based on the minimum water table found over management district 4. The initial water level was set at 3,069 feet based on the irrigated acreage and initial pumping height to the surface of 26 feet, which is the average across management district 4. The depth of the aquifer and initial total acreage surface area make up the two physical parameters needed to define the cone shaped function used to capture the spatial depletion, Figure 1. As discussed in equation (16) above, the gradient of this function describing depletion has implications to the optimal extraction time path, as it defines the magnitude of the marginal cost of lost irrigated acreage due to a change in groundwater height, or our spatial externality.
To estimate crop yields as a function of rainfall and applied irrigation water, we used Kansas State's Crop Yield Predictor tool, which is parameterized for the Colby, KS area (Klocke et al. 2010). We fit functions for corn, sorghum and winter wheat. By running the tool for the full range of water applied to crops and rainfall, we estimated per-acre yield as a function of rainfall, or rainfall plus irrigation water applied. The assumption used here is that an inch of rainfall is equivalent to an inch of irrigation water applied. We fit cubic functions to the yield data. Visually, the yield functions are found in Figure 3. Corn is assumed grown on irrigated acreage and a rotation of wheat, fallow, sorghum, fallow on dryland acreage (Hansen et al. 2012). Thus, the crop yields from dryland are assumed to be 1/3 of the per-acre sorghum yield at a particular Expectations of rainfall play an important role in defining optimal management. The model presented is flexible to various definitions of rainfall expectations and stochastic processes. For clarity of interpretation and limiting the computational burden, we chose to use three, roughly equally likely, levels of rainfall representing low, medium and high amounts of yearly rainfall 10 . We fit an empirical timehomogenous 11 Markov chain process to binned rainfall amounts to match the low, medium, and high levels to observed rainfall at the Colby, KS gauge. The transition probability matrix and details are found in Table 2. 10 Yearly rainfall is the input to the Crop Yield Predictor. From there, a weather generator is used to create growing season weather to estimate crop yield. 11 This implies a stationary distribution and constant transition matrices. By simulating this process for 250,000 years, we estimated the non-conditional probabilities of each state, which are used for our stochastic case, as well as the average of the process for the deterministic treatment to match the conditions for each scenario for comparison.

Numerical Solutions
Optimal policy functions are found numerically by means of stochastic dynamic programming using MATLAB 2013. The optimal extraction decisions as a function of water level are approximated by first estimating the value function, , by means of value function iteration by solving equation (11) for discrete levels of , replacing the optimized values for and repeating until the functions converge within a specified tolerance (Putterman 1994). Creating the transition probability matrix for each realization of rainfall incorporates the Markov chain yearly rainfall process.
Once the optimal value functions are found, the corresponding policy functions can be recovered for any realization of rainfall and groundwater height. The value function and policy function are displayed in Figure 5 and 6.
The derived optimal decision rules are iterated through time starting at the initial groundwater levels. Realizations of rainfall in each year are generated from the i.i.d or the Markov process fitted above to yearly rainfall to match expectations used for transition probabilities. To evaluate gains from optimal management, a myopic decision is made in each year by using the policy of maximizing equation (6) for each year with respect to groundwater withdrawals. These paths through time for the i.i.d case can be found in Figure 6.
The welfare implications of making optimal extraction decisions are estimated by discounting and summing each year's profit as defined by equation (6) from the initial year to the end of the time horizon, in our case 500 years, which is well after reaching a steady state and where future benefits are essentially discounted to zero. The differences between the two welfare values calculated are the estimated gains from optimal management of the aquifer. The results are presented in Table 3 for different specifications of the problem.

Results
We present results from our model applied to the characteristics of groundwater management district 4. The welfare results can be found in Table 3  The results show that including the loss of the spatial extent of the aquifer leads to a larger estimated welfare losses when compared to other modeling assumptions for welfare impact studies over parts of the Ogallala. To put our work into perspective, our estimates are larger than Lee et al (1981) for the Ogallala in Texas at .3%, Nieswiadomy (1985) for the High plains aquifer in Texas at .28%, and similar to Kim et al. (1989) also for the High plains in Texas (1-3.7%). Kim's paper includes endogenous technological change as a function of groundwater heights and therefore similarly provides backstop technologies in a sense, albeit under certainty.
The welfare gains arise from the difference between the optimal decision rules, and the myopic pumpers' choices, Figure 5. The optimal planner internalizes the intertemporal stock externality, while the myopic extractor does not. For the stochastic rainfall scenarios, the realized rainfall determines the corresponding optimal policy function, one for each state of rainfall. As compared to the myopic decision maker, withdrawals following optimal management are about 16% less in low rainfall, 19% less in medium rainfall, and 25% less in high rainfall, varying through groundwater heights ( Figure 5) showing relatively larger savings in better years. When these optimal policy functions are iterated through time ( Figure 6) meaning the decisions are made optimally based on the current rainfall and groundwater state, the groundwater stock is depleted more slowly under optimal decisions, as expected.
The lower optimal extractions, in equivalent groundwater heights and rainfall amounts compared to the myopic pumper, leads the groundwater levels to settle at a steady state height of 2,991 feet on average as opposed to 2,980 feet under myopic extractions for a difference of 11 feet in height at the steady state and between 0 and 17.5 feet through the time-path. With our spatial cone model, where the irrigated acreage is a function of the groundwater height, these steady state groundwater levels translate to an optimal steady state with 118,360 irrigated acres or 31% of the initial irrigated acreage, as opposed to 93,517 irrigated acres under myopic extraction or 25% of the initial irrigated acreage. The sources of welfare losses are the higher pumping costs with deeper water tables and the loss of irrigated acres. The magnitude of these losses is less obvious as they depend on the rate of depletion, the discount rate, and the relative benefits of irrigated farmland as compared to dryland farming.
When comparing welfare across rainfall scenarios, modeling uncertainty increases estimated gains from management, resulting in an increase in welfare gains of .09-.13% over the deterministic welfare gains. The concavity of the one year benefit function implies risk aversion in the sense that irrigation water limits the range of available water in future years, which is preferred over a wider range of possible states as described in Provencher & Burt (1993) and under risk aversion in Knapp & Olson (1996). However, the yield of the backstop technology of dryland farming is less variable (flatter in Figure 3) than irrigated crop yields. In a sense, the transition to dryland farming attenuates the loss of the buffer value of the groundwater, as the crops are more resilient to a range of rainfall then non-irrigated corn would be. Overall, the increase in welfare gains when including uncertainty reflects the additional benefit gained when taking into account groundwater as a steady resource to buffer variation in weather in optimal management.
The groundwater height remains higher through time when taking into account uncertainty in precipitation as the optimal policy consists of slightly smaller extractions on average. The average difference in height is small and only around 1 foot, well within the margin of error of the simulation runs. Clearly, the year-to-year extractions are different in the stochastic cases from the deterministic case, where optimal policy functions imply changing policy rules based on yearly rainfall, Figure   5. So, although the extraction paths are similar in average groundwater height, the stochastic case leads to variable extractions and a range of groundwater levels depending on the realization of rainfall ( Figure 6). The addition of a Markov chain process added little additional welfare gains (.04%) when optimal policy was matched to a stochastic process as described in Table 2.

Discussion
Unlike other approaches, our spatial depletion model captures the gradual loss of access to a stabilizing resource as a result of the overall extraction of the aquifer. The  Table   1).

That modeling the loss of irrigated acreage does not result in significantly larger
welfare losses than what we found was unexpected, although our welfare estimates are larger than those found over similar regions in previous literature. We include a backstop technology, dryland farming, which reduces the gains from management, since in its absence, we would have assumed infinitely expensive alternatives (Koundouri & Christou 2006). Without a backstop, the welfare gains from management are larger at 5.2% in the stochastic rainfall scenario (Appendix Table 1).
We chose to model the dryland crop as a rotation of sorghum, wheat, and fallow in equal proportions. Dryland practices require fallowing some areas and rotations of other crops to maintain soil conditions. We assumed no loss of total farmland, when some land may have to be fallowed permanently after losing access due to its soil characteristics, which may cause us to overstate the value of the backstop technology.
Advances in dryland techniques and drought resistant crops could mitigate some of the negative consequences to losing irrigation access, but even as dryland techniques gain efficiency, so too will the returns to irrigated agriculture and it is unclear if the relative value difference will change much (Steward et al. 2013).
Second, our model assumes that 17% of the initial aquifer surface area is irrigated farmland (Table 1). This is based on the observed number of acres farmed over this district. The total surface area of the usable aquifer, or initial irrigation intensity, may vary widely across aquifers or regions within an aquifer. The results are particularly sensitive to the intensity of farming, or the amount of farming relative to the size of the resource, with a higher percent of acres irrigated leading to larger welfare gains from management (Appendix Figure 1). This suggests that intensely farmed areas or aquifers may have a magnitude of difference in gains from optimal management.
In exploring the implication of the above assumptions, we can say something about when management is likely to have a greater welfare gain. When the value of the backstop technology is small compared to the value of the use with access to the water, the welfare gains are larger. Koundouri and Christou (2006)  We also held prices constant through the time horizon, but adding demand growth for the irrigated crop would certainly increase the welfare gains, similar to what is shown in Brill & Bruness (1994). There may be important price dynamics in large aquifers that support large monocultures, as irrigation water dries up the prices of irrigated crops should increase and the opposite market forces would be in play for dryland crops. These forces would increase the gains from management as the marginal benefits from irrigated agriculture would increase compared to dryland agriculture.
Our motivation in including uncertainty was that groundwater has a value as a buffer to smooth out returns over time (Tsur & Graham-Tomasi 1991;Knapp & Olson 1995).
With a backstop that does not enjoy the stable source of groundwater, we surmised the two aspects together could lead to additional welfare gains from moving away from myopic pumping. Despite the small additional percentage gains from management, a few points can be made from interpreting the optimal rules under uncertainty. The optimal extractions vary year to year, so any management designs to approximate the optimal extractions would need to allow for this variability as rules based on a deterministic optimal solution would limit welfare beneficial pumping in dry years and similarly be too high to induce optimal saving in better years ( Figure 5).
The groundwater stock paths are similar under uncertainty and certainty on average, but vary in the stochastic cases, shown Figure 6. The gains from management could be small and similar to the deterministic case or larger than the average depending on the realization of rainfall, only one of which the resource would actually progress through in reality.
We included a Markov process to investigate if the type of stochastic rainfall process mattered to welfare gains to management, as this is not something that has been discussed to our knowledge. We find that changing the process mattered little to the estimated gains from management, although it changed the optimal policy rules. In our case it slightly increased extractions in each rainfall state, in each groundwater level although the differences are well within the modeling uncertainty. Given that adding uncertainty in general added less than a quarter percent to welfare gains, finer distinctions, such as defining a different process, did not materially matter.
In our model, we simplified droughts with a Markov chain rainfall process where the lengths of the droughts were defined by the transition probabilities. Drought refers to extended periods of lower than normal soil moisture as a result of climatic variables including rainfall (used in this paper), temperature, and even wind speed. Drought measures such as the Palmer Drought Severity Index include a duration aspect, where the progressions of the climatic variables matter and not just the current conditions.
We do not carry over soil conditions or other stock variables, except the groundwater height, that could be affected by the series of yearly rainfall events. This could be done by modeling additional processes in the benefit function and is left for future work.
The addition of the lagged process mattered little to the welfare gains from management and only slightly to the optimal policy functions in our application.
The discussion of rainfall processes leads to an important point about the assumptions leading to the optimal paths. We matched the rainfall processes to the optimal policy rules in each case assuming perfect information about the process that generates the annual rainfall. Any policy rules would be such that the process is not known with certainty. The method we chose to use to estimate optimal management assumes that the decision makers not only act optimally in each period given the realization of rainfall, but also have accurate expectations of the process that generates rainfall. In reality, the processes that generate weather and longer-term climate trends are complex and the distribution of rainfall is not stationary. As such, the extent to which the policy maker's expectation of rainfall coincides with the actual process may have implications to welfare gains, with mistakes or mis-specifications causing possibly large welfare losses.

Conclusion
We introduce a dynamic spatial depletion model of groundwater extraction that incorporates stochastic rainfall and a gradual spatial stock externality, leaving more farmland at the mercy of variable rainfall as groundwater levels fall and less available for irrigation across a smaller area. By building a novel model flexible to rainfall expectations and various stochastic processes, we showed the extent of the importance to welfare and optimal management when including stochastic elements. We also explore the implication of droughts and their persistence on optimal management as it compares to myopic extraction behavior. We find that the addition of randomness and persistence of rainfall does not materially affect welfare gains, largely due to the relatively good yields from dryland farming. Incorporating the spatial depletion of the aquifer added to welfare estimates of moving from open access to optimal management of the resource.

FIGURE 1 Spatial Depletion Model
NOTEu is the radius of the irrigated acreage, x the groundwater height and p is the pumping height, which is the difference between the surface and the height of the groundwater. Θ is fixed given the parameters x, p, and u. Low, medium and high rainfall amounts are shown as vertical lines. The Y-axis represents yield in terms of bushels per-acre as a function of rainfall for sorghum and winter wheat, or rainfall plus irrigation water for corn.

FIGURE 4 Optimal Value Functions
NOTE-This figure plots the optimal value functions for the stochastic (one for each realization of rainfall). The value function represents the net present value of the groundwater resource assuming optimal future management.

Introduction
Natural gas is increasingly being used for electric energy power generation across the United States. This trend is driven by increased domestic production from shale gas as well as regional environmental initiatives.

Background
Natural gas in New England is used by the residential, commercial, industrial and electric sectors. Total yearly consumption has increased from 710,000 to 868,000 MMcf from 2001 to 2014, with the largest percent and total volume gains from electric energy generation (Figure 1). This increased consumption came partly as a result of decreasing prices from a plentiful supply of domestic natural gas due to the technological advances in extraction using hydraulic fracturing and horizontal well drilling (DOE 2013;DOE 2015). For power generation, this has meant the addition of natural gas fired turbines to replace retiring coal and nuclear facilities. Figure 2 shows electric generation by fuel type for New England, showing natural gas replacing coal over the time frame. The total energy generation has remained mostly constant.
The demand for natural gas is highly seasonal and closely correlated with heating and cooling needs. The overall quantity consumed across all sectors peaks in the winter months. However, the electric sector follows a different seasonal pattern. The residential, commercial and industrial sectors use natural gas in higher volumes in the winter for heating, while the electric sector uses more in the summer for meeting cooling demand and less in the winter due to higher spot prices and fuel oil as an available substitute for many generators (Figure 1).
Pipeline and liquefied natural gas (LNG) infrastructure must be sized to meet the needs of the peak demand in the winter, and therefore runs below capacity for the other parts of the year. New England receives natural gas through pipelines that carry it from producing regions in Canada and states to the west such as Tennessee and Pennsylvania. The major pipelines are the Algonquin Gas Transmission, which connects to supplies coming from the South, the Tennessee Gas Pipeline, which runs from the Gulf of Mexico through the producing regions of Pennsylvania and Ohio, and the Maritimes and Northeast pipeline from Canada. In addition, imports of liquefied natural gas (LNG) from overseas come through facilities located around Boston. As there are no major storage sites in New England, the area relies on storage and outflows from the states to the west and LNG imports to meet unexpected demand.
Pipeline customers can be divided into two types: those with firm delivery contracts for a specified capacity, and those that have interruptible contracts for pipeline capacity and the associated gas volumes. Residential and many commercial customers are serviced by local distribution companies (LDCs), utilities who have long-term contracts with the pipeline for capacity to guarantee supply to their customers (DOE 2014). To a lesser extent, industrial users also enter into longer-term agreements for supply, resulting in higher average but less volatile prices. Electric energy generators, on the other hand, generally do not enter into long-term contracts. Their supply is interruptible for a number of regulatory as well as economic reasons explained below.
Therefore, the price the generators pay for natural gas varies greatly and a stable supply is at risk in times of high constraint in the system. The prices across sectors can be seen in Figure 3.
Pipeline owners operate in a highly regulated market, where the price charged for capacity is capped by the Federal Energy Regulatory Committee (FERC) in order to counteract the monopoly power and protect customers (EIA 2015a). However, firm capacity contract holders can sell their owned capacity and volume of gas in an unregulated secondary market. These transactions are made between buyer and seller and are not cleared through a market entity. As laid out in  and , capacity constraints lead to higher spot prices in constrained demand regions and higher price spreads between pipeline nodes.
Therefore, additional rents are made available for firm capacity contract holders with the existence of the secondary market.
In New England, these firm contract holders are the utilities and large industrial users.
Why electric generators have not sought more firm supplies is unclear. Ideas put forward include the availability of alternative fuels for generation (oil), the inability to store gas or electric energy, the ability to pass the higher spot costs along in liberalized electric energy market, the lower average price received from more interruptible contracts in the summer, and issues of counterparty risk due to the heterogeneity in sizes of gas producers and electric energy consumers (Morris 2013;ISO New England 2015). Generators may enter more long-term agreements over time, but it is currently clear that the electric energy sector pays significantly higher prices for natural gas in times of constraint due to the market structure described above.
Constraints are prominent in the winter months when heating needs and natural gas for There is a vast literature of modeling approaches to estimate natural gas price relationships across a multitude of markets. Early approaches focused on end users' costs and substitutability between energy sources using panel and system of equation approaches to identify the elasticities of interest. A review of this early work is presented in Al-Sahlawi (1989). In that time, demand was inelastic over the short term, yet elastic over the longer run, and closely related to the price of substitutes, namely oil. Given the simultaneity and endogeneity issues in estimating price elasticities in a supply and demand system, simultaneous equation models (SEM) and variations of instrumental variable approaches were used for identification of price elasticity and fuel substitution effects (Zellner & Palm 1974).
While these methods lead to a deeper understanding of the long-term drivers of natural gas prices, fundamental factors of supply and demand, such as weather and storage, which can be important at a national and certainly a regional level were largely attributed to unexplained noise. Mu (2007) and ending up with a stationary error term (Murray 1994). For example, this is the case for the Henry Hub natural gas price and the WTI price for crude oil found in Brown & Yücel (2008). Because of this they were able to estimate a long run relationship between the two.
The Henry Hub price is important for many financial contracts and represents the spot price of natural gas in the most liquid market in the US. The difference between this spot price and the price paid by end users in each region has been the topic of a number of studies. Since the 1970s, the market for natural gas has undergone a number of rounds of deregulation. By 1989, wellhead price ceilings had been removed. Resultant regional natural gas spot markets emerged and an interstate network of pipelines developed. In 1992, pipelines were required to provide "openaccess" transportation, separating the pipeline owners' ability to sell the commodity and the transportation capacity together (Apergis et al. 2015). De Vany & Walls (1993) argued that a measure of the extent of integration is to test the level of cointegration between wellhead prices and spot markets separated by a pipeline system. They found that in 1987, before the deregulation, cointegration relationships did not exist between most wellheads and spot prices. By 1991, more than 65% had become cointegrated as a result of the deregulation according to the authors.
Subsequent studies have used this idea and tested the effect of deregulation and connectedness on the co-integration of a number of price series (Walls 1993;Apergis et al. 2015;Arano & Velikova 2010). They each find that deregulation led to a more integrated market where prices are functionally related to each other in the long run.
To what extent pipeline capacity works to increase or decrease the level of integration is not addressed. These papers also do not address a marginal change in capacity on a region's prices or consumption. Arano & Velikova (2010), for example, ask if citygate (the price the local utilities pay) and residential prices are cointegrated within a state.
They find that they are related in the long run post-deregulation supporting the thesis that deregulation has led to better market integration. However, Marmer et al. (2007) use similar cointegration tests to identify constrained regions, whose price time series deviate from long-term relationships with the Henry Hub spot price. Surprisingly, they find that California is the most isolated and that there were no bottlenecks into the Northeast. The constrained conditions have been more recent, while their study's data While including pipeline capacity directly has, to my knowledge, not been included in econometric estimates of New England's regional natural gas prices or quantities, it is often cited as a major driver in elevated prices for the region in recent winter months This paper will build on previous work modeling natural gas prices by including relevant fundamental factors of weather, storage, and alternative fuels while adding and testing the addition of an indicator of pipeline capacity as an important explanatory variable for a regional natural gas market. A description of the data is followed by the econometric methods and results.

Data
The data for the study spans January 2002 to December 2014. Monthly data on natural gas prices and consumption of each state by sector (electric, residential, commercial, industrial) are from the EIA's Natural Gas Monthly reports and obtained through the online data portal. Oil and LNG prices were also obtained through the EIA.
There are multiple prices for natural gas, depending on the market of interest. The city-gate price (Algonquin Citygate for New England) is the wholesale price of gas received by local gas utilities from a pipeline operating company. The sector-specific prices and consumption data comes from surveys reported from individual purchasers and the prices reflect the total cost of the natural gas: the commodity as well as delivery to the end-user (EIA 2015b). These sector prices reflect the type of contractual arrangement that the buyers have with the pipeline or utility for residential, commercial and industrial prices. So, although the Algonquin spot price is a general indicator of the monthly cost of natural gas in the region, it does not reflect the actual price paid by end users.
Pipeline capacity data by state and region is reported annually by the FERC (Federal Energy Regulatory Commission), which is tasked with regulating pipeline use and construction. FERC also approves future projects and release a projection of capacity based on announced or under construction additions and subtractions. Pipeline capacity is represented in MMcf/d, million cubic feet per day. An estimate of pipeline utilization for a month is the difference between total consumption for the month and the maximum possible based on pipeline capacity. For example, for the year 2014 the northeast states had a net import capacity of 120 bcf/month (billion cubic feet) and a consumption of 98 bcf in January, but only 62 bcf in July. For non-firm contract holders, the capacity remaining after firm demand is met reflects their time-varying supply.
Total monthly consumption by sector is used to separate the firm contract demand (local utilities for residential, commercial and industrial uses) from electric sector gas consumption, which occurs under interruptible contracts. The difference between pipeline capacity and firm demand determines the capacity leftover for the electric sector. This is the capacity variation I use to identify the effect of a marginal increase in capacity on the gas market for the electric energy generators. In essence, this would assume price inelastic short-term demand from firm contract sectors. Exact monthly utilization data reported from each pipeline operator and not gleaned from EIA consumption data would be ideal, but it could not be obtained in a consistent manner and is not made available through public reporting.
Previous studies have found strong relationships between the storage, or deviations from normal storage, and the price of natural gas (Mu 2007;Nick & Thoenes 2013).
Natural gas is stored in underground reserves in depleted oil and gas fields and in places where the right geologic properties for storage naturally occur such as salt caves and aquifers. The amount of gas in storage less the base gas needed to keep the storage site operational is known as the site's working gas reserves. The EIA reports the quantity of gas stored by consuming region. This study uses working gas for the eastern consuming region. The extent to which storage affects regional prices within the context of other constraints will be tested in the SEM econometric specification.
Degree days are defined as the difference between the daily average temperature and 65°F. The monthly measure is the sum of the degree days for the month. Both heating and cooling degree days are expressed as positive values. The idea is that cooling demand is determined by how far above the temperature is from 65°F and the opposite for heating. Degree days and differences from their long-run (1981-2010) averages, or degree day anomalies, are strong explanatory variables for natural gas demand and prices. This data comes from NOAA's National Weather Service Degree Day Statistics database, which provides population-weighted heating and cooling degree days by month, state, and region (NOAA 2015).

Model
The economic system can be summarized by a system of simultaneous equations.
There is one equation for each endogenous variable and one identity. Assessing the system in the structural form helps make clear the assumptions of the econometric specifications. (1) (3) The variables are defined as follows: Table 1 Equation (1) is the natural gas demand equation for the electric energy generating sector and equation (2) is the supply equation. Equation (3) is the firm contract natural gas quantity consumed. Equation (4) represents working gas storage. Equation (5) is the identity making explicit the estimate of capacity left for the electric sector being the difference between firm contract consumption and total pipeline capacity.
Including , the measure of capacity for the electric energy sector in the econometric specifications, will test its importance in explaining price and quantity in this market.
The simultaneity between prices, quantities and working gas storage makes identifying the structural parameters ( for exogenous and for endogenous variables) difficult.
For example, to estimate a change in supply as a result of a change in price for natural gas ( ) using OLS, one would have to ensure that the variable for the price ( was uncorrelated with the error term ( in that equation. Given that the quantity and prices are simultaneous, the correlation is guaranteed. To see this, substitute from equation (2) into equation (1) and solve for . This algebra involves dividing equation (1), including the error term, through by a term including the parameter , after the substitution. Therefore, the error term is correlated with consumption, violating an assumption of OLS leading to bias in the coefficient estimate (Murray 2005 (1974). The primary interest of this study, the effect of a capacity expansion on equilibrium price and quantity for natural gas for the electric sector can be addressed with a reduced from equations of the system. Therefore, I estimated the following reduced form equations: The capacity left for the electric sector, , is exogenous to these two equations assuming and do not belong in the equation for , firm contract demand.
Essentially, this assumes firm demand and pipeline capacity are exogenous to the natural gas price and quantity for the electric sector. As the equations are stated in the structural form, there may be an indirect path from , to through working gas reserves, possibly challenging this assumption. Put another way, if changes in working gas reserves affect both firm and interruptible natural gas price supply and demand, which both in turn affect changes in working gas reserves, the firm demand and leftover capacity may be endogenous to equation (6) and (7) above. We will put this path of effects aside for investigation in the estimation of the SEM model, which can deal with this endogeneity.
In order to obtain unbiased and consistent estimates of the reduced form equations above, a number of additional econometric issues need to be overcome. To meet the assumptions needed of OLS to be unbiased and consistent, the error term must be stationary, not correlated with its lags, and normally distributed. Stationarity refers to the property of mean, variance, and autocorrelation of a time series that is consistent across time. Price and quantity time series are notoriously non-stationary.
Regressions between non-stationary variables may lead to spurious results if the two variables happen to follow a similar random path across the observed time period (even if the error term of a regression appears stationary). If this is the case, large values of one series do correlate with large values of the others and vice-versa.
Running OLS would then assume, incorrectly, that the mean is constant which leads to estimating relationships that are meaningless. Regressions between stationary and non-stationary variables also lead to estimation issues, as the non-stationarity is then passed on to the error term. Therefore, I ran tests for stationarity in the form of an augmented Dickey Fuller test, which tests the presence of a unit root, or a coefficient of 1 on a regression of a variable with the lag of itself. These are presented in Table 2.
Those variables where we fail to reject a unit root, implying non-stationarity, are differenced and then re-tested for difference stationarity. The failure to reject a unit root for oil prices and LNG prices implies these are non-stationary series.
The price variables and the measure of leftover capacity are included in natural log form. Therefore, I am estimating a log-log model between those variables and a loglinear relationship between the variables left in levels. Both of these specifications assume a non-linear underlying relationship.  find this sort of nonlinear relationship between the capacity constraint and the price differential between pipeline nodes. How well this type of relationship describes the price and quantity data in the Northeast will depend on the behavior of the error term in the regression.
Namely, after variable transformations, the residuals from the regression should be normally distributed. Using only non-logged data led to a non-normally distributed error term and a fitted model that underestimated the price increases in constrained months. The quantile plots of the residuals of both specifications are found in Figure 4 showing the residuals to be closer to normally distributed (the straight lines) using the log transformations.
In  The results show a significant relationship between the pipeline capacity left for the electric sector after taking into account firm demand and the price of gas for this sector. The coefficient (-.48) means that an increase in capacity left of 1% lead to a .48% decrease in price. The coefficient on degree day anomalies is significant and positive. Degree day deviations remain in level form, so a 1 degree change in degree day anomalies for a month leads to .06% increase in price.

Reduce Form Results
The coefficients on heating and cooling degree days are insignificant. These are highly correlated with the measure of leftover capacity, which is itself a function of degree days through firm demand, equation (3). Including both raises issues of collinearity.
The coefficient correlation matrix reveals that the estimates of the coefficients of HDD, and CDD are in fact correlated with the leftover capacity, our variable of interest (correlations of .82 and .24 respectively). Specifications with combinations of including temperature and capacity variables together and separately are in Table 4 as robustness checks. One would expect collinearity to inflate standard errors of coefficient on both variables, inflate t-values and thus I would fail to reject that coefficients are different then zero more often. The significance of the coefficient on leftover capacity in the model including both leftover capacity and temperature measures remains strong. This means that it is to a larger extent the congestion that drives the price electric generators pay. It also supports modeling the leftover capacity once the firm demand is accounted for when estimating natural gas prices for the electric sector.
That degree day anomalies were significant and not degree days itself may be partly a result of including a seasonal lag, which takes up much of the variation in monthly average temperatures. It is the deviations from normal conditions that provide explanation above and beyond the average monthly conditions. I left the working gas reserve variable out of the regression due to the simultaneity issues. Thus, the estimates of the effect of the temperature variables could be working through changes in working gas as well as through other excluded endogenous variables in the reduced form estimates. In order to separate those effects, one would have to estimate the whole system of equations with identifying restrictions. Table 5 shows the results for the reduced form equilibrium quantity estimates. The coefficients are of the expected sign. The seasonal pattern of natural gas consumption for electricity generation means more is consumed in the summer when prices are low and demand is high and less is consumed in the winter with constrained capacity and high prices (Figure 1). The positive and significant coefficient on CDD reflects the summer cooling demand. Since the equation is in reduced form, the negative sign on degree day anomalies may reflect the effect of the resultant increase in price that go along with those conditions. Again, a system of equations with identifying restrictions would be needed to disentangle these two effects. Lastly, our variable of interest, leftover capacity is significant and positively related to the amount consumed by electric energy generation.

Increase in Capacity Scenario Analysis
The estimates of equilibrium price and quantity allows for estimating the effect of an The regression results for estimating are in Table 6. I use the same process as the previous regressions. First, I fit the OLS model with the covariates, determine the error structure of the residuals, and then re-estimate with GLS accounting for the error structure (ARMA(2,1)). The model is a good fit with a small standard error (.12), which is a 1.1% percent error. The difference between this estimate of firm demand and the pipeline capacity, equation (5), creates the capacity left for the electric sector.

SEM Model
The reduced form models are useful in estimating equilibrium price and quantity, but have their drawbacks. I was unable to include working gas reserves, a possibly important explanatory variable, directly due to endogeneity. The reduced form coefficients are the total, direct plus indirect, effects between the variables in the system. Thus, I may be attributing an effect to degree days, for instance, that is the sum of degree day's direct effect on demand as well as its indirect effect through changes in working gas reserves, if working gas reserves effect supply or demand directly. Additionally, the structural parameters describe the shape of the demand and supply curve. Of particular interest are the price elasticities. Disentangling these effects by estimating the structural parameters that describe the supply and demand curve for natural gas in this market requires estimating the SEM model. To identify the structural parameters, there needs to be included at least as many exogenous variables as there are equations for endogenous variables in the system. I use the temperature series, oil and LNG prices as well as the lags of each endogenous variable as instruments in the 3SLS estimation of the system. Therefore, my model is over-identified with 13 exogenous or predetermined variables, including the lags, to 4 structural equations. The lagged endogenous variables are considered predetermined, therefore not correlated with the current period's error term and exogenous to the current period (Fox 2002;Henningsen et al. 2007). In addition, a minimum of one exogenous variable must be left out of each equation for identification. This is the SEM version of the exclusion restriction in instrumental variable estimates. 12 The 3SLS results are presented in Table 9. Each column presents the results of the identify the path of the effect and the shape of the supply and demand curve. The supply price elasticity is small, negative and significant, which is against what theory would imply. However, very low and sometimes negative supply price elasticities for natural gas, especially over the short run, have been found before, but are not common (Arora 2014). The elasticity results show that both supply and demand are price inelastic over this period. The inelastic demand and supply price elasticities means that small excess supplies or demands are only cleared with a large change in price, which helps explain the volatility of the price for natural gas, especially for interruptible contract holders.
Previous work has found changes in natural gas storage conditions to be significant in describing price movements in reduced form models (Brown & Yücel 2008;Mu 2007), but by including them, each assumes that their measure of storage levels is exogenous to the current price of natural gas. If this is not the case, the coefficients from those regression results are biased. The inclusion of working gas reserves to the SEM as an endogenous variable led to unexpected results. Changes in working gas reserves was insignificant in explaining shifts in the supply curve but positively and significantly related to demand. The results of the working gas structural equation are also unexpected. Changes in reserves are positively related to electric energy consumption of natural gas, which would naturally be a draw on reserves. Working gas reserves are also, curiously, significantly positively related with both HDD and CDD. The effect of degree day anomalies are negative and significant which may reflect the difference from average storage conditions in times of particularly high demand. The results raise more questions then answers for the effect of natural gas storage in the market for the electric generating sector.
The way the variation in working gas explains the dynamics of this system appears to be against what theory would imply. One would expect larger working gas reserves to decrease price by increasing available supply. While I included demands and therefore withdrawals on the reserves, explanatory variables for the variation in additions to working gas storage from producers were not included and could be a direction of future inquiry. Since the seasonal pattern is handled by the seasonal lagged term and the temperature variables, variation in additions to storage above and beyond normal conditions is assumed to be random noise. The gas is being stored to meet expectations of future demand, so if explanatory variables of storage additions are correlated with already included variables to the system, the results could be biased.
Between the curious supply price elasticities and the results of the working gas variable, I do not put much weight on the SEM results at this point other then to reinforce that leftover capacity is a possibly significant driver of supply for this natural gas consuming sector and that the system appears to be price inelastic. These low price elasticities are understandable as supply is fixed with pipeline capacity and demand the result of running large fixed capital.

Discussion
The results of the reduced form regressions show that capacity constraints are a driver of the increase in the price of natural gas in cold winter months. The SEM results support the case for including the capacity constraints as an explanatory variable of the supply of gas to a constrained region such as New England.
Capacity additions represent a shift in the supply curve for natural gas, as identified through the SEM results. The total effect (direct and indirect) of a 1% increase in capacity of prices from the SEM results was similar, but slightly smaller, to that found with the reduced form model, -.44 and -48% respectively. Analytically, these marginal effects on equilibrium prices should be identical assuming the structural model is the correct model of the system. The reduced form model is just a simplification of the structural model. However, the reduced form models are estimated with GLS using maximum likelihood methods, while the 3SLS estimate with GMM. While both are consistent, meaning they converge to the true value of a parameter they no not have identical small sample properties. Additionally, the SEM model explicitly includes the working gas reserves, which could provide a back-door rout for electric sector demand to effect firm demand and thus challenge an assumption needed for the reduced form models to be unbiased. This rout, working gas reserves effecting firm demand, is insignificant in the SEM equations. The two approaches produced similar equilibrium results, but the SEM model allowed for the identification of the shape of the supply and demand curves.
To a large extent, the entire system is driven by temperature dependent demand, so identifying the effect of the capacity constraint above and beyond the effect of the temperature variation was always going to be a tall task. That the variable proposed in this study to represent the leftover capacity for the non-firm contract holders remained significant through many specifications, including ones with other correlated covariates, supports the inclusion of this fundamental factor.
The inclusion of the capacity measure and its significance to explaining both price and quantity provides a path to estimate the effect of increasing the capacity of pipelines Given this non-linearity, pipeline capacity affects the range of prices, reducing the winter spikes to a larger extent in abnormal conditions. Thus, pipeline capacity reduces the price risk of energy supply for the region. However, excess capacity to address winter constrained conditions will come as an additional cost in other times, where the large fixed cost of the system are spread out over the smaller gas volumes that are flowing. A weakness of the methods used in this study is that they do not address directly the possible price changes due to the effect of larger fixed costs of pipeline operation. To what extent the firm demand or interruptible customers share the burden of the increased fixed cost is unclear.
The objective of this paper was to estimate changes in equilibrium prices and quantities of natural gas consumed by the electric sector as the results of a change in pipeline capacity. It is the equilibrium prices for natural gas that describe the per-unit cost of the input to electric generation. These costs are passed along to electric ratepayers. Thus, to estimate the effect of increased capacity on ratepayers, and inform the debate about how the cost of electric energy might change, a reduced form model is well suited. Translating the price for gas that natural gas generators pay to the cost of electric energy to residents and businesses is left for future work. The SEM results clear up the causal mechanism, a shift in supply, and confirm a similar total effect as from the reduced form models. In addition, the SEM model shows that both the supply and demand curves are price inelastic.
This study provides one piece of the puzzle in determining the effects of proposed expansions, but not the whole picture. In the debate about the impact of proposed pipelines, this study offers a way to quantify the effect on the electric energy generators input price and quantity for natural gas. The benefits to the electric energy sector are potentially large by allowing larger quantities of gas to be consumed at lower price for generation. These benefits should be weighed against the costs of the planned routs, the environmental consequences, natural gas price in other sectors and other indirect benefits and costs.     Table 3. The measure of capacity for the electric energy sector remains significant through alternative specifications of the model with combinations of correlated temperature variables 51.371 *** (df = 7; 147) * p<.10 ** p<.05 *** p<0.01 Note: Each column represents a different specification of the model explaining natural gas quantity consumed by the electric sector. The last 2 columns are from GLS regressions with ARMA(2,1) errors specified based on the residuals of the OLS regression. Observation lengths differ due to the inclusion of a seasonal lag. Oil and LNG prices are in first differences in the last column due to non-stationarity issues. 523.714 *** (df = 6; 149) * p<.10 ** p<.05 *** p<0.01 Note: Models for the firm contract quantity consumed are presented in each column. The GLS regression specifies a ARMA(2,0) errors based on minimizing AIC from alternate lag specifications. The results are used to estimate the leftover natural gas capacity for the electric sector. The simulation uses these fitted firm quantity estimates to create the capacity left for the electric sector.  Table 3. The winter scenarios present a range of severities of winters. Price is highest in January. The effect of increased capacity is largest in more severe winters.

FIGURE 1 Natural Gas Consumption by Sector New England
Note: Data is from the EIA. Residential, commercial and industrial uses increase in the winter due to heating requirements, while electric energy uses peak in the summer due to cooling demand.

FIGURE 2 Electric Energy Generation by Fuel Source New England
Note: Data is from ISO-NE and represents the total GWh from each source for the year. Natural gas generators have replaced retiring coal plants. Total consumption has remained mostly constant.

FIGURE 3 Natural Gas Price by Consuming Sector New England
Note: The price data is from surveys conducted by the EIA and accessed through the EIA's data portal. Each sector's price is the per unit price of natural gas including delivery. Residential, industrial and commercial prices follow a different path based on the structure of the local utilities. The electric sector operates under interruptible contracts, which leads to increases in times of pipeline constraint.

FIGURE 4 Log Transformation Q-Q Plot Test Preferred model
Level Model

Log Transformed Model
Note-This figure shows the quintile plot of the residuals from the preferred price model, column 6 in Table 3, when specifications with the dependent and independent variables are in levels compared to in logs. The temperature variables remain in levels in both models. The log transformations create errors that are closer to normally distributed, the strait lines in the plots.

Derivation of Equation 1
The area of a circle representing the irrigated agriculture acreage above the aquifer.
With a fixed initial aquifer area and maximum depth, we can relate geometry of the aquifer to the radius of the irrigated acreage Using these two, we get irrigated area as a function of the groundwater height and geometry of the aquifer.
With the initial conditional known, we can solve for initial radius and : With these we can come to equation (1) in the paper by making at a percent of the initial aquifer area:

Discussion of Equation 8
Here we want to illuminate the model dynamics of including a cone and uniform spatial depletion of the groundwater resource. The speed of depletion is important in the model and is determined by the irrigated area and thus a function of the height of groundwater; therefore we compare the speed of height change in our model to the speed of height change in a bathtub model.
Proposition 1: given the same demand for water and a decreasing path of groundwater height. The change in groundwater height is larger in the bathtub model than in the spatial depletion model.
We start with defining per acre water demand constant to give a fair comparison between models where is the per acre water demand in the bathtub and spatial depletion model at a given height defined by the price of water P. The total amount of water extracted is a function of the height of the water as is the change in groundwater height. The total water extracted decreases in the spatial model because there are less irrigated acres as the height of groundwater decreases. This is realized through the spatial depletion function or equation (1) in the main text. The rate at which groundwater height changes, the equation of flow, is also a function of . refers to the spatial depletion model and refers to the bathtub model. and Equation (4a) describes how the percent of irrigated area of the of the spatial depletion model changes. Increase in height cause an increase in irrigable area, and a decrease in height causes a decrease in irrigable area at this rate described in in (4a). Taking a difference in equation (2a) and (3a) we find Equations (6a) and (7a) state that as the aquifer is depleted the change in height is less in the spatial the depletion model at the rate given in equation (4a) and because the spatial component becomes smaller and is in the denominator of equation (2a) the recharge acts as a buffer to loss in height in the spatial depletion model. As the depletion parameter, , gets closer to one the difference converges to zero. As decreases so does , which drives the differential to be positive.

APPENDIX FIGURE 1 Welfare Gains to Management as a function of farming intensity
Stochastic Scenario NOTE-Welfare gains from management are sensitive to the assumption of farming intensity which defines demand for groundwater and its marginal value. As farming intensity is increased so to are the estimated gains to management.

ARIMA Reduced Form Model
The I in ARIMAX stands for the level of integration, or how many times the dependent variable must be differenced to be stationary. Using the stationary series, having first differenced the non-stationary variables, an ARMAX of the following form is estimated.
Where is the lagged, period, dependent variable with coefficients . Seasonal lags enter in the same manner, but lagged 12 periods on our model. These terms handle the autoregressive, AR p , process.
are coefficients on the moving average, MA q , process on pervious residuals of a lag length . X is a matrix of exogenous regressors with coefficient estimates . The appropriate length of the lag is a choice of the modeler, but can be guided by evaluating the autocorrelation correlation function (ACF) and partial autocorrelation function (PACF). From there it is common to use measures of AIC and BIC can determine the most appropriate lag lengths.
I estimate the reduced form equations for both price and quantity based on the criteria described above, which turned out to be the lag structures found in the OLS residuals and specified in the GLS regressions from Table 3 and 5. The price equation is of the form ARIMAX(1,0,0)(1,0,0) 12 and the quantity equation of the form ARIMAX(2,0,0)(1,0,0) 12 , where ARIMAX(p,d,q)(P,D,Q) m . The second set of parentheses specifies the seasonal terms. M is the seasonal period lag. The results are presented in Table 8.
The results are similar to that found in the GLS regressions with an increase in leftover capacity effecting prices negatively. The coefficient is larger in magnitude, implying a 1% change in capacity leading to a .6% change in price. Degree day anomalies are significant and of similar magnitude. The quantity equation similarly shows a larger effect of the leftover capacity on the equilibrium quantity consumed.
The results confirm the results of the GLS regressions. A benefit of the ARIMAX models is that predictions can be made that take into account the time series of the error structure, as opposed the GLS results, which use the unbiased estimates of the coefficients to estimate values, but does not take into account the recent error term innovations in making step ahead forecasts. Therefore, once values for the exogenous regressors are provided, future values can be assessed under both current and increased capacity. Similarly to the policy simulations above, first I estimate the firm demand, using an ARIMAX model, calculate the remaining capacity and estimate the price and quantity for electric energy generation equations. I estimate 2 years ahead, with the time series of the severe winter scenario. The results are in Figure 6. p ** p *** p<0.01 Note: Reduced from equation of the natural gas price and quantity for the electric energy generating sector are presented. The lag structure is based on minimizing AIC with the inclusion of the exogenous covariates and a 12-month seasonal lag.

APPENDIX FIGURE 2 Natural Gas Prices for Electric Energy Generation Equilibrium Price and Quantity ARIMAX Models
Price Quantity Note: Price and Quantity projections from the ARIMAX model in Table 8. Months count ahead from January, 2015 to December 2017. The red line represents an increase of 450 MMcf/d of pipeline capacity added at the stat of 2016, or 12 months ahead in the figure. The increased capacity leads to lower prices and higher quantity consumed.