Containing the Risk of Phosphorus Pollution in Agricultural Watersheds

Phosphorus (P) is an essential nutrient to boost crop yields, but P runoff can cause nutrient over-enrichment in agricultural watersheds and can lead to irreversible effects on aquatic ecosystems and their biodiversity. Lake Erie is one prominent example as this watershed has experienced multiple episodes of harmful algal blooms over the last decades. Annual P loads crucially depend on yearly weather variations, which can create the risk of years with high runoff and excessive nutrient loads. Here we apply stochastic modeling to derive sustainable management strategies that balance crop yield optimization with environmental protection, while accounting for weather variability as well as weather trends as a result of climate change. We demonstrate that ignoring annual weather variations results in mitigation efforts for environmental pollution that are largely insufficient. Accounting explicitly for future variations in precipitation allows us to control the risk of emissions exceeding the P target loads. When realistic risk targets are imposed, we find that a package of additional measures is required to avoid P over-enrichment in the Lake Erie watershed. This package consists of a substantial reduction of P inputs (approximately 30% for different accepted risk levels), adoption of cover crops throughout the nearand mid-century, and cultivation of less nutrient-intensive crops (30% more soy at the expense of corn). Although climate change reinforces these conclusions, we find that the accepted risk level of exceeding P target loads is the predominant factor in defining a sustainable nutrient management policy.


Introduction
Agricultural and industrial development has led to nutrient over-enrichment or eutrophication in surface waters in the last decades [1]. This anthropogenically induced abundance of nutrients can favor cyanobacteria, and result in harmful algal blooms (HABs). HABs, in turn, negatively affect the habitat of many fish species and can result in excessive fish mortality through oxygen depletion. Other ecosystem services where HABs have negative outcomes include the use of fresh water for irrigation or drinking water, and recreational activities. HABs are a threat to the integrity of water bodies all over the globe, including Lake Victoria in Africa [2], Lake Erie in the US and Canada [3], Lake Taihu in China [4], the Baltic Sea in Europe [5], or the Caspian Sea in West Asia [6].
Lake Erie is the shallowest among the Great Lakes of North America and is particularly vulnerable to eutrophication and resulting HABs. In this study, we concentrate on the Western Lake Erie Basin (WLEB), which spreads across three states and spans nearly seven million acres. We focus in particular on the Maumee River watershed, which is the largest watershed of Lake Erie. This is a highly agricultural area (≥80%) and is shown to be the largest contributor of agricultural P to Lake Erie. Several studies have confirmed that the Maumee River is the dominant source of P loading causing cyanobacteria outbreaks in the WLEB, and bloom sizes over the past decades have been closely related to external P loading from the Maumee River [7][8][9]. Over the last decades, continuous efforts were made to protect the biotic integrity of Lake Erie [10]. Point source policies had a large effect on soluble P concentrations in the 1970s and 1980s [11]. While agricultural policies have historically been centered around sediments, the scope of agricultural policies has broadened starting from the 1980s [12]. In order to combat the observed re-eutrophication since the mid-1990s, the Objectives and Targets Task Team for Lake Erie recommended a P load reduction of 40% in the Western Basin tributaries with respect to the 2008 water year loads, based on a range of modeling load-response analyses [13].
In addition to the intensification of agriculture, climate change poses serious threats to aquatic ecosystems [14]. Rising temperatures have a non-linear effect on P nutrient uptake [15], and can favor cyanobacterial growth and occurrence of HABs [1]. Climate change will also concur with hydrologic changes through modified patterns, intensity, and duration of precipitation events. The growing intensity of precipitation will increase both nutrient runoff from agricultural fields and groundwater discharge, ultimately leading to more nutrient enrichment of the surface waters. In fact, climate studies have predicted the trend of increasing P loads under climate change [14,16], and the increasing frequency and magnitude of rain events have been indicated as the most important factor in the increase of soluble P loads [12,17]. As P transfers will increase, the implications for agricultural practices to mitigate these future emissions are still unclear. The P target loads chosen by the Objectives and Targets Task Team are all based on the assumption that precipitation patterns do not change in the future [13], which is debatable given the anticipated changes in the Lake Erie ecosystem due to climate change [18]. Concerns that climate change will exacerbate the frequency and magnitude of extreme weather events, therefore, necessitate the development of adaptive nutrient management strategies.
Several simulation studies exist that evaluate the effectiveness of nutrient management in reducing nutrient loading under various climate scenarios [19,20]. Nutrient management models are typically deterministic and have limited ability to account for weather variability in the definition of optimal farming practices, balancing farmers' income with environmental impacts. Although deterministic models are able to present a diversity of nutrient management practices under different scenarios, this diversity can only be obtained under conditions of perfect foresight. In real-life conditions where the amount of precipitation in future years is unknown, the ability of these models to define optimal nutrient management strategies is more limited. This study contributes to the literature by evaluating how weather variations in a changing climate, unknown at the moment that decisions are taken, affect sustainable nutrient management practices. To this purpose, we draw on methods from stochastic optimization, which have recently been applied in the field of agricultural hydrology and irrigation management [21,22]. A small number of stochastic optimization models exist that deal with nutrient pollution, for instance for Lake Balaton in Central Europe and Erhai Lake in China, but these studies do not consider climate change [23,24]. In this work, we develop a decision-making framework that provides optimal packages of sustainable agricultural practices according to different preferences regarding accepted risk levels. Specifically, we integrate a risk-averse strategy that accounts a priori for weather variability and allows to set reliability targets for environmental protection. In addition, we include a deterministic, risk-neutral strategy that provides optimal farming practices based on expected precipitation levels, indifferent to the variability in precipitation that inevitably will affect the P loading into the surface waters. As random precipitation patterns play a crucial role in nutrient runoff, this framework allows us to quantify unplanned excess P loads that result from ignoring weather variability. Furthermore, the developed framework enables us to benchmark the package of nutrient management practices that results from a model that ignores risks versus a model that anticipates risks. Finally, as climate change will alter P loadings in the agricultural watershed, we quantify its effect on future phosphorus transfers, and analyse to what extent the package of nutrient management practices has to change to contain the risk of excessive phosphorus pollution.

Materials and Methods
In order to define a sustainable nutrient management policy, we construct a modeling framework around three building blocks: (i) an economic component that captures profits from agriculture, (ii) a soil component that describes the phosphorus cycle, and (iii) an environmental component that reflects the decision-maker's perception of environmental risk. The framework provides a set of optimal agricultural management practices, i.e., land allocation for different crops, fertilizer (P) application to the soil, and the use of cover crops in the off-season. Decisions on land allocation, fertilizer application, and cover crop adoption all contribute to the nutrient dynamics in the soil, and ultimately to the nutrient emissions into the surface waters. The framework takes the perspective of a regional decision-maker who protects both the interests of farmers and the environment. The decision-making framework maximizes the profits of farmers, but provides different strategies for environmental protection. In the risk-neutral strategy, the regional decision-maker uses the best available estimates of the emission rates to make sure that the expected P loads do not exceed the P target loads. In the risk-averse strategy, the regional decision-maker accounts for the inherent uncertainty in emission rates due to random precipitation events, and sets reliability targets for not overshooting the emission limits. The regional planner using the risk-averse strategy aims for high reliability of the environmental protection, or equivalently, low accepted risk levels of emissions beyond the P target loads.

Economic Component
The profit π(t) of crop production on a representative hectare can be decomposed into contributions π i (t) of the considered crops i, i.e., corn, soy, and wheat, each under 4 different levels of conservation tillage. The profit function can be written as where δ i (t) is the land allocation indicating the share of land per hectare allocated to each of the crops. Crop profit π i is a function of three decision vectors: fertilizer application F(t) in kg/ha, cover crop adoption θ(t) representing the share of a hectare used for cover crops in the off-season, and land allocation δ(t). The argument ω is a realization of a random vector capturing the random annual amount of heavy precipitation that results in random emissions of P into the surface waters. The profit per crop can be written as The profit function accounts for the income from crop production and multiple cost terms. The income from crop production depends on the crop yield Y i , the price p i , and the variable costs C v,i . The variable costs are modified by the level of conservation tillage indicated by the crop residue level R i . Prices and variable costs are decreasing over time at a decreasing rate (see Table A2 in Appendix B). The input costs are divided into two parts, those affected by tillage practice (second line in (2)) and those that are not (third line in (2)). Among costs that are influenced by tillage practices, in addition to the variable costs C v,i which vary with yield, we account for fixed costs C f , chemical costs C ch , and labor costs, which are modeled as the product of labor price p l and labor hours H i . Since conservation tillage reduces tilling operations, the parameters a 1,i through a 4,i are positive parameters that are calibrated to indicate the cost-effectiveness of conservation tillage. Note that while conservation tillage is generally associated with lower input costs, chemical costs are higher compared to conventional tillage [25,26], so the sign before a 3,i is positive. The costs unaffected by tillage practices include the cost for fertilizer application with P price p f , the cost of animal manure M with application cost p m , and the cost of cover crops with annual cost per hectare denoted by C c . Per hectare crop yields (in kg/ha) are defined as Crop yields are a function of soluble P reserves in the soil P ss and the crop residue level R i . Note that we do not consider annual yield growth related to improving farming practices.

Soil Component and the Phosphorus Cycle
Common nutrient management practices that can reduce P emissions from agricultural land include precision nutrient management, conservation tillage, and cover crops. Conservation tillage is defined as any practice that leaves at least 30% crop residue on soil surface [27]. It has been widely promoted as an effective way to combat P runoff by trapping P in soils on farm fields. Despite its widespread application in the US, the benefits of conservation tillage are debated with regard to its effects on crop yields and efficiency in controlling P runoff [26,28,29]. Conservation tillage likely reduces soil erosion, and hence reduces sediment-attached P loss from agricultural fields [30]. There is mounting evidence, however, showing that conservation tillage increases the export of total soluble phosphorus [28,31]. Since soluble P is almost 100% bioavailable [32,33], the effect of controlling P export by conservation tillage could be offset by soluble P loss [34,35]. Planting cover crops has also been shown to be an effective strategy to reduce P emissions from agricultural land [36][37][38]. By absorbing the residual P, cover crops will reduce the number of nutrients potentially available to be transported to surface waters. Cover crops can also reduce attached P since they reduce soil erosion.
The dynamics of P in the soil system are modeled via two state-equations for sediment attached P and soluble P. The stochastic state equation for sediment attached P, P sa (t, ω), is given by The state variable P sa is assumed to be uniformly distributed in the representative hectare of the farm. The carry-over coefficient ζ 1 ∈ [0, 1] represents the proportion of existing attached P reserves carried over from the last period to the current period. The parameter τ A is the proportion of manure that goes into the attached P pool. Conversion between sediment attached P and soluble P is captured by the final two terms in equation (4). The fraction of attached P that converts to soluble P each year is given as φ A , and the fraction of soluble P that is converted to attached P each year is given as φ S . Finally, the sediment attached unit area P load E sa (t) from fields into river systems can be written as The parameter γ sa i (t, ω) is the uncertain emission rate at time t with ω capturing the uncertainty of nature, in this case related to the amount of annual heavy precipitation. The emission rate describes the share of P stock that results in the annual P load. Conservation tillage traps sediment attached P by reducing soil erosion, which is captured by the term 1] and b A accounting for the effectiveness of conservation tillage in reducing sediment-attached P loss. The effectiveness of cover crops in reducing attached P loss is represented by the term (1 − c A ).
The dynamics of soluble P in the soil, P ss (t, ω), are represented by the following stochastic state equation The carry-over coefficient ζ 2 ∈ [0, 1] represents the proportion of existing soluble P reserves carried over to the next year. Fertilizer application F i adds to the stock of soluble P in the soil, and τ S is the fraction of annual manure application that turns into soluble P. Annual consumption of P by crops is given by (4), the conversion between sediment attached P and soluble P is governed by the parameters φ A and φ S . Finally, the unit area P load of soluble P, E ss (t), into river systems is given by (7) which depends on the random emission rate γ ss i (t, ω). Conservation tillage in this case increases soluble P runoff, with b S the effect of conservation tillage on soluble P emission. Cover crops act as a mechanism to reduce soluble P emissions, whose effectiveness is captured by the parameter c S .

Risk-Neutral Nutrient Management
The risk-neutral regional decision-maker aims to maximize the profits of farmers while ensuring that the P loads do not exceed the emission limits. In order to do so, a risk-neutral decision-maker makes use of the best available estimate of the emission rates for sediment attached P and soluble P. This problem can be formulated as the maximization of expected profits under dynamic constraints (4) and (6), and P loading constraints that ensure that P loads do not exceed the imposed P target loads η sa and η ss . The optimization problem can thus be written as s.t.
(4) and (6) In reality, emission rates can diverge from their expected values, but the risk-neutral decision-maker does only account for nutrient emissions under expected precipitation conditions.

Risk-Averse Nutrient Management
In every year with major algal outbreaks in Lake Erie, big spikes in discharge from the Maumee watershed were observed [39]. A major part of the uncertainty in P loading can therefore be attributed to stochastic weather events, which are captured in the model by adjusting the nutrient emission rates. A risk-averse regional decision-maker aims to hedge against the risk of environmental costs associated with heavy precipitation events. A possible approach is to define nutrient management policies under worst emission conditions, but this strategy is likely to cause high losses in agricultural profits. Since meeting the emission targets can be infeasible under certain circumstances, or can come at an exceedingly high cost, the regional decision-maker can also impose a target reliability level, i.e., the probability with which P emissions must be below the P target loads. In this case, the optimization problem can be relaxed to a chance constrained problem with reliability level ρ sa and ρ ss , with 0 ≤ ρ sa , ρ ss ≤ 1, and can be written as s.t. (4) and (6) The risk-averse decision-maker maximizes expected profits and requires that the emission targets are met with a high probability. This problem formulation is numerically intractable, yet equivalent to a two-stage optimization problem, where strategic actions (land allocation, fertilizer application, and cover crop adoption) are complemented with adaptive actions in case the actual P loads exceed the P target loads. These adaptive actions can correspond to a penalty for excess emissions, or to the cost for water treatment [40,41]. For more details on how to solve the chance-constrained optimization problem, we refer to Appendix A.

Uncertainty under Different Climate Scenarios
Weather plays a crucial role in P loads predominantly through rainfall. In fact, the weather has been estimated to account for about 50% of the higher soluble P loadings in the period 1996-2011 [12]. In the soil component of the proposed framework, variable emission rates are the considered source of uncertainty and are contingent on the amount of yearly heavy precipitation. The P stock in the soil depends on the sequence of emission rates, and therefore on the sequence of yearly precipitation. To compute the probability of meeting the emission constraints, we need to construct uncertainty scenarios over the considered time period. In other words, we create realizations of the emission rates γ sa i (t, ω s ) and γ ss i (t, ω s ), where ω s represents a specific realization of the underlying uncertainty, and where 1 < t < N t and 1 < s < N s .
Climate change may dramatically alter model parameters like crop growth rates, nutrient emission rates, and the persistence of HABs [1,42,43]. In this work, we are mainly interested in how climate change will affect precipitation patterns and the corresponding emission rates. We consider three climate scenarios: a stationary climate, an intermediate scenario corresponding to RCP 4.5, and a worst-case scenario corresponding to RCP 8.5.
To generate uncertainty scenarios in a stationary climate scenario, we sample from the empirical probability distributions of the emission rates of sediment attached P and soluble P. The distributions are derived from annual phosphorus loading data for the Maumee River for water years 1975-2015 and historical soil test levels [44,45], and are depicted in Figure 1. Note that the empirical distributions are asymmetric with longer right tails. From the time series data, we observed a very high correlation between the emission rates of sediment attached P and soluble P, and we assume here full correlation between both emission rates. To generate uncertainty scenarios under the RCP 4.5 and RCP 8.5 climate scenarios, we make use of the Maumee River discharge projections in the 21st century [17]. Using the SWAT hydrological model and an ensemble of global climate models, the river discharge has been estimated under both climate scenarios for the near-century (NC: 2010-2039), mid-century (MC: 2040-2069), and late century (LC: 2070-2099). We assume here that the relative increase in annual discharge can be used as a proxy for the relative increase in the emission rates. The annual discharge is projected to increase by 6.5% (12%) in the near-century, and 7.7% (14%) by mid-century under RCP 4.5 (RCP 8.5). We model the non-stationarity of the emission rates by shifting their respective distribution functions according to the corresponding relative increase, both in the near-and mid-century.

Results
The decision-making framework described in the previous subsection is used over the time window 2016 to 2050, under three climate scenarios and four different preferences for environmental risk, i.e., the risk-neutral strategy and the risk-averse strategy with reliability levels equal to 85%, 90%, and 95%, respectively.
The initial values for the P stock state variables are set equal to their historical observed value. Target loads for both sediment attached P and soluble P are set according to the recommendations by the Objectives and Targets Task Team. Specifically, their recommendation is a 40% reduction of P loads with respect to the 2008 water year loads, aiming to achieve a bloom size no greater than that observed in 2004 or 2012, 90% of the time [13]. We study the three major crops in the watershed region (corn, soybean, and wheat), and assign four conservation tillage categories within each crop. Specifically, we assign four values for the residue levels R i , according to minimum (0%), low (35%), high (75%) and maximum (100%) [46]. Crop prices and variable costs are assumed to follow their historical trends after 2015. Cropland proportions are also assumed to be bounded by a ±50% deviation following historical patterns for the Maumee River basin, using county-level harvested crop acre data (USDA-NASS, 2016). A list of the model variables and parameters can be found in Tables A1 and A2.

Emissions for Different Reliability Targets and Climate Scenarios
The likelihoods of meeting the target loads for sediment attached P and soluble P are very different (Figure 2). Using the risk-neutral strategy for sediment attached P, the reliability of meeting the P target load is above 90%, and the risk-averse strategy further increases the reliability up to approximately 100% (Figure 2a). The recommended emission limit for sediment attached P is undemanding, but this is however not the case for soluble P. In the remainder of the paper, our focus is therefore on soluble P loads. For this type of P loading, the risk-neutral strategy only achieves a reliability level of 52% as can be observed from the cumulative distribution function of soluble P loads under a stationary climate (Figure 2b). By introducing reliability targets in the risk-averse strategy, reliability levels up to 95% can be achieved, which are observable by the strong shift of the distribution function to lower P load values. Evaluating the magnitude of excess P loads, we find that the risk-neutral strategy results in excesses between 43% and 48% on average above the P target load under RCP 8.5 and the stationary climate scenario, respectively ( Figure 3). In other words, disregarding the stochastic nature of the phosphorus cycle would produce emissions above the 2008 water year load levels approximately every second year. Conversely, acknowledging the stochastic nature of the phosphorus cycle allows for setting reliability targets in the riskaverse strategy, which reduces the magnitude of the average excess P loads substantially. Specifically, with a reliability target of 95% the expected excess P load is marginally above the P target load (<5%), and this result is robust under the three climate scenarios.
The regional planner with a risk-neutral approach towards environmental risk determines solution pathways for P input, cover crop adoption, and land allocation, using the best available estimate of the emission rates. Using the risk-neutral strategy, the P target loads for soluble P are only achieved with a probability of 52%, which is largely insufficient to fulfill the objective of a 40% reduction with respect to the 2008 P loads [13]. The risk-neutral strategy actually generates solutions that result in expected P loads beyond the 2008 water year loads biennially. As a consequence, using the solutions as specified by the risk-neutral strategy in a realistic environment with random weather variations provides an unacceptable level of environmental protection. Magnitude of excess P loads. When P loading occurs beyond the P target load, these excess emissions are on average between 43% and 48% higher than the P target load in the risk-neutral strategy (i.e. these events happen each year with a likelihood of 48%). When a reliability target of 95% is imposed, the excess P loads are only marginally above the P target load.

Tradeoff between Farmer Profits and Reliability of Environmental Protection
The results from the former section call for the introduction of reliability targets and for models that account explicitly for weather variability in the definition of reliable solution pathways. In line with the loading targets of the Objectives and Targets Task Team, we evaluate reliability levels between 85% and 95%. Meeting the reliability targets for P loading comes however at a price in terms of farmer profit loss (Figure 4). With respect to the risk-neutral strategy, the risk-averse strategy results in a 15% profit loss for an 85% reliability level in a stationary climate, and a 28% profit loss for a 95% reliability level under RCP 8.5. Farmers will endure substantial profit losses to meet the emission targets recommended by the Objectives and Targets Task Team, and climate change aggravates this effect moderately. To make sure that farmers adopt the reliable solution pathways recommended by the risk-averse nutrient management scheme, farmer losses will need to be compensated through crop price adjustments, subsidies, or compensations by other interest groups. Private farmer profits are represented on the vertical axis relative to the profit per hectare of the risk-neutral strategy in a stationary climate. Note that private farmer profits do not include the environmental and social costs due to intensive farming. We consider different risk preferences on the horizontal axis, and the three considered climate scenarios are indicated by different colors. Meeting reliability targets for P loading comes at a high cost. A reliability level of 95% corresponds with a profit loss of 23% under a stationary climate, and a loss of 28% under RCP 8.5.

P inputs, Land Allocation, and Cover Crops
Considering P inputs aggregated over the period 2016-2050, the risk-neutral strategy under a stationary climate recommends on average P inputs of 22.5 kg/ha/y. The fertilizer application rate in the Essex region at the shore of Lake Erie was 26 kg/ha/y in 2006 [47]. In other words, the risk-neutral strategy prescribes in this region an average reduction of P inputs by 16% with respect to 2006. Application rates vary, however, and the USDA WLEB CEAP assessment for 2012 reported an overall average application rate of approximately 19 kg/ha/y. Following the risk-neutral strategy, P inputs should drop by up to 7% with respect to a stationary climate when we account for climate change (see Figure 5a under the risk-neutral strategy). However, when environmental risk is considered in the risk-averse strategy, P inputs need to be reduced by approximately 30% with respect to the risk-neutral strategy under a stationary climate. The effect of introducing reliability targets dominates the effects of climate change. Interestingly, with a reliability level of 95%, the recommended P input can be higher under RCP 8.5 than under a stationary climate, but this result needs to be considered in conjunction with the adoption of cover crops, which is considerably higher under RCP 8.5 (Figure 5c). However, it stands out that the solution pathways produced by a deterministic perspective on the phosphorus cycle in the risk-neutral approach provide largely insufficient mitigation efforts to control nutrient pollution.
Evaluating the dynamics of P input over time provides more insight into the underlying processes (see Figure 5b for P input recommendations under the RCP 8.5 climate scenario). During the initial years, all strategies recommend low P input levels to consume the very high levels of P stock in the soil. From the 1960s to the 1990s, the recommended level of P inputs for optimal crop nutrition was markedly higher than the amount that crops could remove [48]. This trend was exacerbated under the promotion of conservation practices in the WLEB during the past decades, especially in the form of sediment-attached P, creating the legacy P problem [49][50][51]. Due to the high levels of P stock, optimal nutrient management involves drawing down this legacy through the use of P already stored in the agricultural systems. At the end of the time period, we observe again low levels of P inputs. This could be partially related to the numerical effects of terminal conditions over the considered time window. However, in the risk-averse strategy with a reliability level of 95%, the decrease of P inputs by mid-century is clearly observable. The higher the reliability level, the more conservative P inputs become by mid-century, even though in combination with higher P inputs during the near-century.
The risk-neutral strategy recommends the adoption of cover crops only in the first years to help reduce the P stock in the soil, and this recommendation holds for the different considered climate scenarios (Figure 5c). Conversely, the risk-averse strategy recommends cover crop adoption throughout the entire time period. In a stationary climate, cover crops are initially used in 100% of the agricultural land, and cover crop adoption declines to approximately 20% by mid-century. Under RCP 8.5, higher levels of cover crops are maintained throughout the near-century, and cover crop adoption drops to approximately 30% by mid-century.
The risk-neutral and risk-averse strategies allocate around 2% of land to wheat, both in the near-and mid-century and under different climatic assumptions (Figure 5d). However, we observe that corn is substituted by soy by mid-century, and this effect is more prominent under RCP 8.5. The risk-averse strategy allocates in general more land to soy than to corn as part of the strategy to control nutrient emissions and move away from the most nutrient-intensive crop, corn. Regarding conservation tillage, all strategies recommend high residue levels throughout the near-and mid-century. (c) (d) Figure 5. Levers for pollution control. (a) The total P input over the considered time period should be reduced by approximately 30% for the different considered reliability levels with respect to the risk-neutral strategy. The risk-neutral strategy in a stationary climate requires 22.5 kg/ha/y on average. (b) P input trajectories are shown under RCP 8.5. P input profiles are smoothed using a 6-point moving-average to reveal better the P input trends. All profiles of P input start with lower levels to consume the existing level of P in the soil. When a reliability of 95% is imposed, the P input is higher in near-century (NC), while the P input drops substantially by mid-century (MC) to adjust for the increased amount of precipitation. (c) Under the risk-neutral strategy, cover crops are only used initially to get rid of the high P concentrations in the soil. The 95% reliability strategy prescribes however the continued use of cover crops up to 2050. Under RCP 8.5, high levels of cover crops are maintained throughout the near-century. (d) Crop allocations exhibit interesting dynamics between NC and MC. In the risk-neutral strategy, corn is substituted by soy and the effect is substantially more prominent under RCP 8.5. Using the 95% reliability strategy, less corn is produced in favour of soy.

P Stocks for Different Reliability Targets and Climate Scenarios
As mentioned above, the nutrient management for sediment attached P is less critical than for soluble P. The solution pathways obtained from the risk-neutral and risk-averse strategies all result in gradually decreasing stocks of sediment attached P, with a reduction of 25% by 2050 for the risk-neutral strategy and approximately 30% reduction for the risk-averse strategies (Figure 6a-c). Apart from minor differences, these results hold for the three considered climate scenarios. Soluble P in the soil follows a very different trajectory and sustains a major drop over the first years, mainly thanks to reduced P inputs and also through continued use of cover crops in the risk-averse strategies. In order to meet the 40% reduction of P loads, P stocks in the soil need to be lowered by 55% with respect to 2015 conditions according to the risk-neutral strategy. However, the risk-averse strategies prescribe a reduction of P stock in the soil of 70% by mid-century (Figure 6b-d). (d) Figure 6. P stock trajectories. (a-c) Sediment attached P stock in a stationary climate (a) and under RCP 8.5 (c). Both trajectories show only minor differences, and the stock of sediment attached P continues to decrease until 2050. (b-d) Soluble P stock in a stationary climate (b) and in RCP 8.5 (d). After an initial steep decline, the soluble P stock stabilizes around a value depending on the imposed reliability target. In case of the 95% reliability strategy, there is a clear difference between near-century and mid-century P stock values.

Conclusions and Discussion
The objective to reduce P loads in Lake Erie by 40% with respect to the 2008 emission levels requires a coherent package of nutrient management policies. By means of a new modeling framework that computes P loads originating from a representative hectare, we demonstrate that the mitigation efforts to avoid excessive nutrient pollution are largely insufficient if weather variability is not accounted for. A risk-neutral approach, ignorant about deviations from the expected precipitation, results in nutrient management recommendations that overshoot approximately every second year the emission limit by almost 50% on average. As variations in annual precipitation have a sizeable effect on P loads, nutrient management schemes need to account for these uncertainties in order to reliably meet the emission targets. Introducing reliability targets makes a substantial impact on nutrient management. Specifically, compared to the risk-neutral strategy, a robust finding is that nutrient management cognisant of weather uncertainties requires significantly lower levels of P inputs (∼30%), higher adoption of cover crops throughout the near-and mid-century, and land allocation for less nutrient-intensive crops. Climate change further mediates these outcomes, and the effect of increased precipitation is mainly evident by mid-century.
Meeting the emission targets with high levels of reliability comes at a price. Nutrient management strategies that reliably limit P loads imply profit losses for farmers up to 28% under RCP 8.5. However, the tradeoff between private farmer profits and the reliability of environmental protection does not take into consideration the social and environmental costs that result from unsustainable nutrient management. In fact, private farmer profits could also be benchmarked with the environmental cost of exceeding the P target loads (for instance through the cost of clean-up actions), as well as the effects on fisheries and the recreation sector. The chance-constrained optimization problem, which arises from setting risk targets, can be solved by introducing a penalty for excess nutrient emissions. The resulting penalty function problem, also known as a two-stage optimization problem, incorporates the social and environmental costs in the search for optimal agricultural policies. Once the strategic nutrient management is put into practice, emissions can still exceed the safety limits in rare events, and these excess emissions are penalized by means of a fine, tax, or cost for clean-up. In other words, the two-stage stochastic program maximizes net profits, i.e., the private farmer profit adjusted for environmental costs, and internalizes the externalities of intensive farming (social and environmental costs) in the definition of optimal nutrient management strategies.
Proposing nutrient abatement strategies is a difficult task due to the complexity of landscape processes and multiple uncertainties in input data, model structure, and model parameters [14]. In light of these difficulties, the modeling framework developed in this work does not capture the full range of complexities and is an aggregate model, attempting to provide general guidelines. There are several avenues to make the framework more realistic. For instance, the framework can be made geographically explicit, and provide nutrient management strategies considering the spatial variation of parameters in much more detail. Furthermore, a number of farming practices have not been considered in this work, and if added could provide a completer perspective on nutrient management. As an example, subsurface drainage and irrigation can have an important effect on surface runoff and P transfers in agricultural watersheds [52]. Also, the effect of nutrient prices could be included, as they have a non-negligible impact on nutrient concentrations in surface waters [3]. Finally, solving stochastic optimization problems is computationally demanding, and a balance needs to be struck between model complexity and computational feasibility to solve this type of problem. Model relaxations or computational strategies need to be explored to compute optimal management strategies using a reasonable time horizon in combination with a sufficiently large set of scenarios.
This study proposes a nutrient management strategy that balances environmental protection with agricultural profits. Although the results show that environmental protection is possible given the necessary mitigation efforts, agricultural profits are estimated to go down up to 28% under a worst-case climate scenario. Since environmental costs are spatially separated from where agricultural profits are made, incentives will be necessary to achieve the required change in farming practices. In this paper, we are not specific about how this behavioral change can be accomplished. Both rewards and penalties can be used to motivate the adoption of the recommended nutrient control strategy. Farmers can be compensated for their profit loss, or instead they can be made responsible for the environmental damage through taxes or fines. A tax or fine that corresponds to the social and environmental cost of excess P loads is also the interpretation of the two-stage optimization problem as defined in Appendix A. Yet, in the context of nutrient pollution it is not straightforward to hold polluters accountable. Nutrient concentration measurements in the watershed do not allow us to allocate responsibility to specific farms. Nonetheless, P stock levels are measurable, and monitoring the difference from recommended P stock levels (see Figure 6) can be a useful tool to strengthen environmental governance.

Appendix A. Reformulation of the Chance Constrained Optimization Problem
In order to solve optimization problems with probabilistic constraints, we typically need to resort to numerical methods. In particular, in this study only the empirical distributions of the emission rates are available, which results in discontinuities in the probabilistic constraints, making the problem numerically intractable. In the following steps, we reformulate (9) in order to attain a tractable form of the problem. Under some technical conditions [53], it can be demonstrated that the chance constrained problem in (9) is asymptotically equivalent with the following penalty function problem (4) and (6) P ss (t, ω) ≥ 0, P sa (t, ω) ≥ 0, P ss (0, ω) = P ss,0 , P sa (0, ω) = P sa,0 , with Q sa (t, δ, F, θ, ω) = max{0, E sa (t, ω) − η sa } and Q ss (t, δ, F, θ, ω) = max{0, E ss (t, ω) − η ss }. The variables Q sa (t, δ, F, θ, ω) and Q ss (t, δ, F, θ, ω) represent the excess P loads beyond the P target loads and are the optimal solution of the second stage problem min y with x ∈ {sa, ss}. The second-stage problem introduces a buffer variable y that needs to compensate for the excess P loads once the random weather events realize. The chance constrained problem defined in (9) is therefore equivalent to a two-stage problem where the management practices δ, F, and θ are anticipative or strategic decisions that are taken before a specific weather pattern is realized. Anticipative actions take into account expected excess P loads E[Q sa ] and E[Q ss ]. In case the P target loads are not met, a recourse action is taken proportional to the excess emissions, i.e., α(Q sa + Q ss ), with α the cost per unit of excess emission. The equivalence between the problem with chance constraints and the two-stage stochastic optimization problem is established through the relationship between the reliability level ρ in the chance constraint problem and the penalty α in the two-stage problem. Specifically, for a chosen reliability level ρ sa and ρ ss , the corresponding value of α can be found numerically.
The problem defined in (A1) contains the non-smooth functions Q sa and Q ss , which can introduce numerical complications. We circumvent this issue by introducing auxiliary non-negative functions µ sa (t, ω) and µ ss (t, ω), which leads to the following problem max δ,F,θ,µ sa ,µ ss (4) and (6) We assume independent realizations over time of the uncertainties γ ss i (t, ω) and γ sa i (t, ω). However, at each time t we assume full correlation between both random variables. Taking the expected value over a stochastic process is typically intractable, and therefore we approximate the expected values by a sample mean, over a set of N s scenarios. The problem can now be formulated as max δ,F,θ,µ sa ,µ ss (4) and (6) µ sa (t, ω s ) ≥ E sa (t, ω s ) − η sa , µ sa ≥ 0 µ ss (t, ω s ) ≥ E ss (t, ω s ) − η ss , µ ss ≥ 0 P ss (t, ω s ) ≥ 0, P sa (t, ω s ) ≥ 0, P ss (0, ω s ) = P ss,0 , P sa (0, ω s ) = P sa,0 . (A4) The chance-constrained problem with discontinuous constraints in (9) is reformulated as a penalty function problem with continuous constraints, which is in its current form numerically tractable.

Appendix B. Notation
The notation used throughout the paper is summarized in Tables A1 and A2.

Notation Description Unit
π(t, δ, F, θ, ω) profit at time t for a representative hectare $/ha π i (t, δ, F, θ, ω) profit at time t corresponding to crop i $/ha δ(t) land allocation vector at time t for all crops i unitless F(t) P input vector at time t for all crops i kg/ha θ(t) cover crop vector at time t for all crops i unitless Y i (P ss , R i ) yield of crop i at time t kg/ha P sa (t, ω) stock of sediment attached P in the soil at time t kg/ha P ss (t, ω) stock of soluble P in the soil at time t kg/ha E sa (t, ω) P load for sediment attached P at time t kg/ha E ss (t, ω) P load for soluble P at time t kg/ha γ sa i (t, ω) emission rate for sediment attached P at time t unitless γ ss i (t, ω) emission rate for soluble P at time t unitless ω uncertainty vector representing random annual precipitation NA Table A2. Model parameters. Parameter values are based on [54] and references therein, and informed by expert assessment. For R i , 4 values corresponding to 4 different conservation tillage practices are given. All other sets with 3 values correspond to the 3 considered crops: corn, soy, and wheat. Both crop prices and variable costs are decreasing with time according to p i (t) = p i (0) · exp(η p,i /∆ p · (1 − exp(−∆ p t))) and C v,i (t) = C v,i (0) · exp(η c,i /∆ c · (1 − exp(−∆ c t))), with η p,i = {−0.008, −0.007, −0.007}, η c,i = {−0.01, −0.009, −0.01}, ∆ p = 0.03, and ∆ c = 0.03.