Probabilistic simulation of phosphorus loss using the Vermont P-index: a bottom-up field to watershed approach

Abstract

Anthropogenic nutrient loading from land use, especially agriculture, is a major threat to waterbodies worldwide. Efforts to govern nutrient pollution are increasingly based on simulation modeling for research, evaluation, and regulation. This study develops a novel approach to simulate nutrient losses from agriculture applied to the Lake Champlain basin in the US state of Vermont. The Vermont Phosphorus-Index—a farm-based empirical model regularly used for site evaluation—is scaled up to the basin level with high-resolution geographic data and probabilistic estimation of unknown parameters and management practices. Results are comparable with analyses using more data and computationally intensive tools. Important insights into basin-wide management include: (1) nutrient-management planning can significantly reduce P losses in a livestock-agriculture-dominated watershed by re-distributing manure applications from areas of high loss to low loss; (2) hotspot identification from geographic data alone may be deeply complicated by high underlying heterogeneity of soil phosphorus; and (3) probabilistic modeling using simple, field-scale models is a potentially useful complement to complex watershed process models. Findings suggest that currently available best-management practices will likely be insufficient to reach reduction targets in the most impaired sub-watersheds. Reductions of agricultural land use and herd size, particularly in intensive dairy operations, may be necessary.

Introduction

Agricultural runoff is the leading cause of nutrient pollution in waterways worldwide (Dupas et al., 2015). Many freshwater bodies are phosphorus-limited, thus nutrient enrichment can profoundly impact their ecology and function (Sterner, 2008). Globally, the scale of phosphorus loss into waterways is several times higher than the geologic background rate, with potentially irreversible impacts on freshwater ecology (Carpenter and Bennett, 2011). The resulting eutrophication also can harm the health, well-being, and economic vitality of communities near affected waterbodies. Driven by excess phosphorus levels, cyanobacteria blooms result in waters with low levels of oxygen, higher turbidity, and dangerous concentrations of cyanotoxins. These conditions can drive die-offs in aquatic animals and losses of native plant species, increase the costs of treating water for human consumption, and threaten the scenic and recreational and scenic values of freshwater lakes (Lewis, Wurtsbaugh, and Paerl, 2011). Exposure to toxins from blue-green algae also may be linked to long-term chronic human health problems (Carmichael, 2001; Torbick et al., 2018) and deaths of bathing pets (Hauser, 2019).

In the state of Vermont, several waterbodies have been declared ‘impaired’ by phosphorus levels, including Lake Champlain, the sixth largest freshwater lake in the USA. While phosphorus runoff is driven by increases in developed land and impervious surfaces, agriculture remains the largest source of phosphorus into the lake, especially in the most impaired sub-watersheds (US EPA, 2016). Water-quality impairment of Lake Champlain has led to numerous negative impacts, including loss of tourism revenues, decreased home prices, and reduction in fish populations (Gourevitch et al., 2021; Voigt, Lees, and Erickson, 2015), as well as less-tangible losses to human well-being (Lake Champlain Management Conference, 2010).

Nutrient runoff from agriculture is considered ‘non-point’ source pollution, occurring at many locations distributed across the landscape. This makes direct measurement of waterbody loading infeasible. Instead of direct measurement, computer modeling plays a central role in assessment and governance (Lane et al., 2006; Wang et al., 2020). While the final levels of pollution or water-quality impairment may be directly measured by water-quality monitoring programs, assigning the responsibility for this pollution to different land uses is accomplished through complex, data-intensive watershed models, especially the soil and water assessment tool (SWAT) (Arnold et al., 2012). These same tools are utilized for assessing different management scenarios and determining what sectors can make the most cost-effective reductions.

In contrast, field-scale and validated models are used to assess on-farm phosphorus loading and management options. This includes phosphorus-indexes, or P-indexes, which are spreadsheet models that can be completed by farmers, extension agents, or consultants. Models such as the P-indexes have become integral parts of both voluntary and mandated nutrient management plans across the USA (Mallarino et al., 2002; Sharpley et al., 2017). In Vermont, the Vermont P-index is a required component of nutrient management plans, which are mandated since 2016 for all farms with more than 50 acres or more than 50 dairy cows (Agency of Agriculture, 2018; USDA NRCS, 2021).

In this paper, we develop a ‘bottom-up’ method for estimating basin-wide phosphorus losses by estimated the P-index for individual farm fields. We link the Vermont P-index with a high-resolution dataset of individual farm fields in a four-county Vermont portion of the Lake Champlain watershed. Being far less computationally intensive than SWAT, it is easier to explore a range of management scenarios and perform calculations on extremely high-resolution land-cover data. Crop fields are analyzed as complete units rather than pixels, closely matching the units on which farm management decisions are made. Privacy limitations preclude access to the state's legally mandated P-index data, so this is a ‘next best’ approach to simulate phosphorus reduction scenarios at farm-scale, aggregate the results through policy and management scenario analysis, and identify potential phosphorus hotspots in impaired sub-watersheds where farm retirement may be the only option to achieve reduction targets.

Background

Modeling agricultural nutrient losses

Computer modeling is a crucial tool for governing the environment, including agricultural runoff. Models are used for many purposes relating to agricultural pollution to waterbodies, including assessing watershed-scale impact, strategic and tactical planning, farm-scale nutrient planning, payments for conservation programs, cap-and-trade design for nutrient management, and even as evidence in lawsuits (Radcliffe, Freer, and Schoumans, 2009). The same model or set of models may be used to diagnose a problem, then design and implement interventions and finally to assess their effectiveness.

Many types of models are used to estimate pollution runoff, ranging from mechanistic models built-up from laws of physics and chemistry to empirical models which utilize coefficients developed from experiments, and mixed methods in between. Models are developed based on different time and spatial scales—from hourly or daily time-steps to annual averages—and from the plot scale to large watersheds. Each approach has distinct strengths and weaknesses, with trade-offs evaluated by the ultimate end use.

For example, dynamic models with large mechanistic components and short time-steps are common for creating regulations and supporting research (Radcliffe, Freer, and Schoumans, 2009). Because these models directly encode a theoretical model, they are useful for advancing scientific understanding. Dynamic models are also generally robust enough to be used for basin-scale environmental assessments. More simplified empirical models, on the other hand, are more user-friendly and have smaller data requirements. The design of these models is also more flexible, where new experimental results can be used to alter key coefficients even when their theoretical basis is unclear. These properties are beneficial for on-farm nutrient management planning (NMP) and assessment of conservation initiatives (Radcliffe, Freer, and Schoumans, 2009; Sharpley et al., 2017).

Watershed-scale nutrient loading is usually investigated utilizing spatially explicit models that are mostly mechanistic. Most popular is the SWAT. First released in 1994, SWAT has become the tool of choice for modeling erosion, streamflow, and nutrient loss from agriculturally dominated watersheds (Gassman et al., 2007), accounting for nearly half of published catchment-scale modeling research (Fu et al., 2019). SWAT is based on a daily time-step and explicitly models the physical processes of plant growth, nutrient uptake, and nutrient transport. By leaving several important parameters to be fit to available watershed data, SWAT can be calibrated to acceptable accuracy for any agriculturally dominated basin.

Among empirical, field-based models, the P-index has become the most popular approach for ‘site assessment’ to help farmers and natural resource professionals determine which farm fields require mitigation measures (Lemunyon and Gilbert, 1993; Sharpley et al., 1994). P-indices are built for utilization in nutrient-management planning using basic information about the physical characteristics and management of farm fields. By utilizing data and calculation tools (such as worksheets and spreadsheets) that farmers and extension agents can easily acquire and use, the P-index method is particularly suited for assessing immediate field-level management actions.

While dynamic models can be calibrated to many different watersheds, the design of P-indices has become increasingly detailed and locally specific. The first P-index was designed for the whole USA, and took an additive approach, using the weighted sum of scores on several different P-loss risk factors. More recent P-indices divide phosphorus loss risk factors into two categories: (1) transport factors, such as erosion and runoff; and (2) source factors, including soil test phosphorus and the quantity and method of recent manure or fertilizer applications. Many recent P-indices also multiply source risks by transport risks (Gburek et al., 2000), providing a much more sensible view of nutrient loss; a field that has no phosphorus source will have no phosphorus loss, no matter how high its transport potential, and vice versa.

Some newer P-indices separate source and transport factors specific to different loss pathways, utilizing regionally specific empirical relationships between management practices, environmental characteristics, and phosphorus loss. For instance, the Vermont P-index multiplies source and transport factors specific to five different pathways: particulate P from soil and manure, and dissolved P from soil, manure, and fertilizer. Because these P-indices are based on experimental and observational data, rather than explicit simulation of the processes governing phosphorus losses, these models are limited to a very small geographic scope and well-researched farming techniques, but can be just as accurate within these limitations as difficult to calibrate, data-intensive watershed models (Sharpley et al., 2017).

An important difference between these two approaches is data needs and thus user preferences. SWAT requires data for many different geophysical features and processes, including elevation, land cover, plant growth patterns, tillage implements, and daily data for rainfall, solar radiation, wind speed, and relative humidity. Additional streamflow and water-quality data are needed for manual or automatic calibration of SWAT (Arnold et al., 2012). The P-index, on the other hand, requires detailed information on management decisions and soil characteristics for individual fields. While a farmer might need some help from an extension agent or technical advisor for determining soil erosion or initial determination of soil characteristics, almost all information needed to run the model is easily accessible to farmers.

Our approach is therefore an inversion of the normal approach. We utilize a tool built for individualized site assessment to examine watershed-scale nutrient-loading. By aggregating simulated results on individual fields up to the watershed scale, we can approach aspects of farm-management that are not available in spatial data as heterogenous and probabilistic, rather than assuming uniformity. For instance, rather than assuming that tillage practices are uniform according to some set of conditions, we can specify tillage and any other farm practices according to a set of conditional probability distributions. This gives a more detailed perspective on uncertainty at the field level, and also may give a more accurate assessment of the average value of P losses; average P loss for a field may not be the same as P loss under its average management parameters. This approach also connects the watershed-scale view of the problem to the farm manager's perspective and direct experience.

Furthermore, our approach has other advantages which may make it useful for analysis of agricultural water-quality policy. For example, a bottom-up approach may be better adapted to handle the considerable uncertainty inherent in these problems. Being far less computationally intensive than SWAT, it is easier to explore a wide range of management scenarios and perform calculations on extremely high-resolution land-cover data through a probabilistic approach. We are then able to explore unseen but potentially highly influential variability in factors such as pattern tile drainage and soil test phosphorus. In addition, the inputs of the P-index are already required elements of legally mandated nutrient management plans, potentially allowing for validation through public records requests. The ad-hoc, empirical nature of the model also makes it simple to adjust its coefficients in response to new data in a constantly evolving landscape of research on water-quality best-management practices (BMPs).

Study site

This study was conducted in the Vermont portion of the Lake Champlain basin in the northeastern USA. Lake Champlain is the sixth largest freshwater body in the USA, draining a basin of 21,326 km² in the states of Vermont and New York, and the Canadian province of Quebec (Fig. 1). Phosphorus-driven eutrophication has been a major environmental concern in the region since the 1970s. Programs targeting point-source pollution from wastewater treatment facilities have been successful, reducing phosphorus pollution by over 80%, but phosphorus concentrations in lake segments drained from Vermont land have remained steady or increased (Smeltzer, Shambaugh, and Stangel, 2012). In 2002, the Vermont Department of Environmental Conservation developed a total maximum daily load (TMDL) plan for reducing phosphorus loading. However, after little state progress and following litigation, the US Environmental Protection Agency (EPA) established newer, more ambitious targets in 2016.

Figure 1. Study area.

The EPA's revised TMDL calls for approximately two-thirds of Vermont's load reductions to come from agriculture, though they estimate agriculture as only producing 41% of phosphorus loading (US EPA, 2016). In response to the litigation, and anticipating the revised TMDL, the Vermont State Legislature passed the Clean Water Act of 2015 (Act 64) Act 64 directed the Vermont Agency of Food Agriculture and Markets to ‘Required Agricultural Practices’ regulation. This regulation came into force in 2018 and aimed at minimum standards for protecting water quality including buffer zones near surface waters. Additionally, all but the smallest farms are required to develop nutrient management plans for all their fields, including using the Vermont P-index and altering management if fields scored too high.

Vermont's agricultural sector is dominated by dairy farming, and the large volumes of manure and imported feed inherent to this sector make reducing phosphorus loading exceptionally difficult (Wironen, Bennett, and Erickson, 2018). Previous studies have investigated the various challenges of meeting TMDL targets for phosphorus, including modeling studies conducted for development of government policy and for more pure research goals (Gaddis and Voinov, 2010; Ghebremichael, Veith, and Watzin, 2010; Meals and Budd, 1998; Medalie, Hirsch, and Archfield, 2012; Mendelsohn, Swanson, and Isaji, 1997; Mendelsohn and Rines, 1995; Seltzer and Wang, 2004; Wironen, Bennett, and Erickson, 2018). Most importantly, the TMDL and reduction allocations for Vermont were developed using the BATHTUB model for within-lake phosphorus dynamics and SWAT for phosphorus losses from land-use, both in the current situation and reduction scenarios. A ‘Scenario Analysis Tool’ was also developed using SWAT (Tetra Tech Inc., 2015) and a literature review to project possible reductions from BMPs in agriculture, forestry, and developed land. These estimates were used to assign load reduction targets to different pollution sources and are utilized in policymaking and assessment by the VT state government.

On the research side, basin-scale studies include Wironen, Bennett, and Erickson (2018), a University of Vermont study that examined the phosphorus mass-balance of the basin as a whole, finding that imports of new phosphorus in feed and fertilizer for dairy farms exceed exports of phosphorus in milk by over 4 to 1. These results suggest that field-scale mitigation measures cannot meet phosphorus reduction targets on their own. Researchers have been particularly interested in the dynamics of the Missisquoi Bay watershed, the most impaired lake segment (US EPA, 2016). Ghebremichael, Veith, and Watzin (2010) utilized SWAT to identify critical source areas in the Rock River sub-watershed, showing that phosphorus loading is highly spatially concentrated. The majority of the load (58%) from the Rock River is generated by fields growing corn (Zea mays L.), despite these fields covering only 17% of watershed area. Winchell et al. (2011), also using SWAT, similarly found that in the Missisquoi watershed, cropland growing corn accounted for only 10% of the land area but about 35% of the phosphorus load. In contrast to Ghbremichael et al., they found substantial hotspot impacts within land-use categories. For example, Winchell et al. found greater than a tenfold increase in phosphorus losses from poorly drained compared to well-drained corn fields, and a fourfold difference between the flattest and steepest corn fields. Similarly, Gaddis and Voinov (2010) built a customized hydrologic model of phosphorus loss in the St. Albans Bay watershed, another significantly impaired lake segment in the northeast corner of the basin. Like other analyses, they estimate that cropland growing corn accounted for a disproportionate quantity of phosphorus loading, with 46% of total phosphorus load estimated to come from only 22% of the land area.

Methods

In this study, we conduct simulations of phosphorus loss from farm fields from four counties within the Lake Champlain basin. These counties account for over 68% of the agricultural land in the Vermont portion of the Lake Champlain basin and over 99% of agricultural land in these four counties are within the basin. Both empirical and simulated variables are used to populate the Vermont P-index, a field-scale model of phosphorus loss used for NMP in Vermont (Jokela, Tilley, and Faulkner, 2020). The Vermont P-index is publicly available as a spreadsheet application, which we converted to a script in the Python programming language (Version 3.8). The Vermont P-index returns a quantitative score with a qualitative interpretation. For instance, a score below 20 is ‘Very Low’ while above 100 is ‘Very High.’ The results can be interpreted as total phosphorus loss, adjusted for availability to cyanobacteria blooms, by dividing by 80 to yield lbs P/acre/year (or dividing by 89.8 to yield kg P/ha/year). Internal calculations of the P-index model can likewise be used to estimate total P losses into waterways.

This modeling approach has several components, outlined in Figure 2. Polygons representing the spatial extent of crop fields in the Lake Champlain basin were developed by the University of Vermont Spatial Analysis Laboratory. These land areas were then linked with elevation models, soil type, and waterway extent from publicly available geodata. Management and other non-observable parameters were modeled probabilistically, detailed in Tables 1 and 2. A description of key modeled parameters is provided below, with more detail in Appendix A.

Figure 2. Tasks performed for this study.

Table 1. Sources for parameters estimated from available data for the P-index model

Notes:

1. Retrieved from https://geodata.vermont.gov/.

2. See: http://www.omafra.gov.on.ca/english/engineer/facts/12-051.htm.

Table 2. Description of simulated parameters for the P-index model

Key modeled parameters

Soil erosion

The Vermont P-index requires a soil erosion rate as a key input. Our model calculates erosion through an implementation of the Universal Soil Loss Equation. The topographic (LS) factor is calculated using raster digital elevation models from the Vermont Center for Geographic information, using flow-routing algorithms from the Whitebox Tools open-source geoprocessing library (Lindsay, 2016), and formulas described by Desmet and Govers (1996). The soil erodibility (K) factor was extracted from the USGS soil map, and a grid of rainfall erosivity (R) factor values was built from data from the EPA Low Erosivity Waiver Erosion Factor Calculator application program interface (US EPA, n.d.). These data were averaged for each delineated crop field to estimate potential erosion by field. This ‘RKLS’ factor for crop fields is held constant across all management scenarios, as it represents ‘potential erosion,’ i.e., soil loss when the field is left completely bare. A lookup table was built for crops, crop sequence, and tillage type to extract crop cover (C) and practice factors (P), which allows erosion to be calculated for each simulated crop sequence and tillage type.

Crop rotations

Expected crop rotations were provided along with the field shape delineations. When we compared these crop rotations with total crop acreage for corn and hay within the study area (NASS, 2021), we found that the provided crop rotations implied a greater ratio of land growing corn to land growing hay in any given year. To correct for this imbalance, corn–hay rotations where corn was represented in more than half of years were altered to add 1 year of hay and subtract one of corn.

Manure applications

Total manure applications were first calculated by county using livestock numbers from the USDA Census of Agriculture (NASS, 2021). Total P applied in manure lbs/year for each county was calculated as:

(1)$${\rm manure\_applied} = \left({n_{{\rm dairy}\,\,{\rm cows}} + \displaystyle{{n_{{\rm other}\,\,{\rm cattle}}} \over 2}} \right)\times 70$$

where n represents the number of each animal. Farmers were assumed to spread manure on hay fields at 80% the regular rate received by corn fields. Each scenario in each county included a series of tuning runs to calibrate the manure application rate such that the total manure P applied matched the total manure P available for that county. Tuning was required for scenarios including NMP, where manure applications are revised if P Index values are too high.

Soil test phosphorus

Values for soil test P were simulated based on data from 6293 soil test samples collected in western Vermont farm fields between 2013 and 2020 and analyzed by the University of Vermont Agricultural and Environmental Test Lab. While zip codes for each mailed sample were available, the data were too sparse to draw conclusions about the spatial distribution of soil test P. Instead, we combined all data and experimented with different distributions. A lognormal distribution was chosen as fitting the underlying data best. See Appendix B for more details.

Soil test P values are extremely influential to P losses at the field scale and are not linked to the spatial data used to run the model, calling for careful treatment in the model. Our process for simulating soil test P was as follows. First, a distribution was simulated from the empirical soil test P data. From this distribution, we drew 10,000 random values, and binned these values into integer-based categories. We then took the mean of each category yielding 84 unique soil test P values. All crop fields, in all scenarios, are simulated with each of these values the same number of times. This allows for analysis of ‘hotspots’ and changes across scenarios to be conducted based on averages for both fields and soil test P results. Average values of results for each crop field are calculated as the average across values of soil test P used, weighted by the number of observations in that bin.

Running the model

For farming practice parameters, several scenarios are defined. These represent different proportions of crop management practices dependent on crops and different field characteristics. For each scenario, an array of values is drawn for each simulated variable equaling the total number of fields. Then for each of 252 iterations of the scenario, these values are randomly shuffled and the values applied to fields in that order. In total, 252 iterations allow each field to be simulated three different times with each of 84 soil test phosphorus values.

A few scenarios were investigated, including: (1) base scenario, representing practices prevalent in 2015; (2) base + NMP scenario, with the same practices but where farmers eliminate or reduce P applications if the calculated P index exceeds certain thresholds; and (3) BMP scenario, representing near-universal, but imperfect uptake of BMPs. All practice scenarios were also run with an alternative set of field shapes where all areas of each field which were within 12.2 m of a watercourse were converted to vegetated buffers. Lastly, all field shapes are simulated under the conditions where they were left fallow and allowed to return to natural vegetation to investigate the impacts of targeted land retirement.

More details on the implementation of the model, including instructions for accessing the source code and description of variables within each scenario, are available in Appendix A.

Results

Impacts across different scenarios

Figure 3 summarizes the simulated phosphorus reductions achieved from different changes in farming practices. The BMP scenario, when combined with nutrient management, exceeds TMDL targets for P reductions in four out of eight lake segments examined, exceeding these targets by a wide margin in Isle La Motte and Shelburne Bay. Extending buffers up to 40 feet away from all surface water yields a 10% reduction in P losses, and together with BMPs and NMP brings the Main Lake and Otter Creek segments into compliance with the targets set for agriculture. In the Missisquoi Bay watershed, less than two-thirds of targeted reductions are achieved in the maximum scenario without land retirement, and only 75% are achieved for Lower Lake A (Table 3).

Figure 3. Total proportion of reductions in P loss from Vermont watersheds under different scenarios.

Table 3. Fractional phosphorus reductions by lake segment from various scenarios

Numbers in bold show where potential changes in agricultural management fail to reach the level set by the TMDL.

These results indicate that achieving targeted reductions in P loss from farmland will almost certainly involve either a radical transformation of farming practices or a substantial shrinking of agriculture, both in terms of animal numbers and land footprint. While an empirical model such as the Vermont P-index is very poorly equipped to simulate novel and radical changes to farming practices, the results of our simulations can help to understand efficient targeting for land retirement.

Of the examined components, NMP has the largest spatial variability. NMP reduces predicted P loading by less than 15% in the Isle La Motte Lake segment, and over 26% in the Missisquoi Bay segment. Higher reductions are seen in watersheds that are partially or totally in Franklin county, which has the highest livestock densities in the state. These reductions primarily come from reductions in manure losses from corn fields with high-soil erosion. Reductions from the suite of BMPs are fairly uniform, ranging from 23% in the Missisquoi Bay watershed to nearly 32% in Isle La Motte. Implementing buffers show modest variation in improvements, with no clear pattern in the variation, as shown in Figure 3.

Mode and quantity of loss

Fields growing corn had substantially higher and more variable P loss than fields in hay in the base scenario. Field-level average P-index values were 77 for corn vs 32 for hay, with 90% intervals of 19–158 and 8–71, respectively (Fig. 4). Losses from manure, however, were somewhat higher for hay (corn P sub-index 16, hay P sub-index of 19). This likely reflects two factors. First, most manure spread on corn fields is incorporated through tillage in the base scenario, while nearly all manure applied to hay fields is surface applied and thus more vulnerable to runoff (despite having far more vegetative cover which can catch that runoff). The incorporated manure P is less-likely to be lost in run-off and more of it ends up in the soil P pool—where it counts as soil loss if lost to subsequent erosion. Additionally, the NMP scenario causes some corn fields to have lower average manure applications than hay fields. Corn fields with medium-to-high potential erosion will sometimes generate a P-index that requires that no manure be applied to the field. Figures 5 and 6 show the relationship between losses from soil and losses from manure for crop fields, with Figure 5 showing results under the baseline scenario, and Figure 6 showing results for NMP and BMPs. Note that NMP allocates manure away from fields highly vulnerable to P loss, creating an inverse relationship between manure P losses and soil P losses in corn fields in the NMP + BMPs scenario.

Figure 4. Histograms of field-average P-index values by crop type.

Figure 5. Relationship between losses from soil and losses from manure for all draws from all crop fields in the base scenario.

Figure 6. Relationship between losses from soil and losses from manure for all draws from all crop fields in the base + NMP scenario.

Spatial vs soil test targeting

Our results suggest that some watersheds, including the Missisquoi Bay watershed, cannot meet phosphorus reduction targets through agricultural management changes alone. Therefore, it may be necessary to reduce the footprint of agriculture or reduce the amount of corn grown within these watersheds. If measures such as land retirement or converting to continuous hay are undertaken, choosing the right fields to target could yield substantial gains in efficiency. Risk of P loss from a field is a combination of innate site characteristics and impacts of past and current management. Innate site characteristics such as soil type, proximity to water, and rainfall erosivity are well-measured and do not need to be simulated within the model. Soil test P, however, is highly right-skewed, apparently log-normally distributed. Detailed spatial data is not available for soil test P.

Figure 7 compares average modeled reductions from field retirement for results binned by location vs results binned by soil test phosphorus for all fields in the sample. The blue distribution represents calculations of mean P-indices conditional on each value of soil test phosphorus, averaged across all fields growing the crop. The red distribution represents calculations of the mean conditional on each field, averaged across all values of soil test phosphorus. Both soil test P data and field geographic data are available, but data linking the two are not.

Figure 7. Distribution of P-index reductions from retiring a field based only on soil test phosphorus or based only on spatial characteristics.

As Figure 7 illustrates, a very small number of soil test P values create extreme hotspots, but the distribution of average P-index scores is highly compressed. This is due to the log-normal distribution of soil test P values (see Appendix 2). Though the absolute spread is wider, the distribution is more unequal for results grouped by field delineation. The Gini coefficient for P-index grouped by fields is 0.48, but only 0.41 for results grouped by soil test P.

Factors determining P-index and P loss of fields

To determine sensitivity of the P-index to various simulated factors, we performed a linear regression with field-level fixed effects on all instances of any field growing corn in the base model run. The results of this regression are given in Table 4. Most of the variation between runs on the same field can be explained by a few simulated variables, including soil test P, cover cropping, tillage type, P added in manure, and whether manure was injected (if the field was no-till). Coefficients for tillage and no-till with manure injected are interpreted relative to a base case of no-till with manure left on the surface.

Table 4. Linear regression coefficients for factors influencing field-level P-index and total phosphorus lost for fields growing corn

Regression contains field-level fixed effects.

Most importantly, doubling soil test P increases modeled P loss by 21%. While this is somewhat inelastic, this also implies that BMPs will not be able to overcome high legacy P, which can result in soil test P levels an order of magnitude higher than those found in a typical field. The model shows substantial P reductions from combined manure-injection/no-till systems.

Recall that the P-index represents total P lost, with erosion losses adjusted downward for P availability in surface waters. All dissolved P and manure particulate P is fully available, while some P in soil particulates is unavailable to algae. Because of this, some variables reduce total P loss more than they reduce the P-index because they primarily reduce P loss through erosion of soil particles more than through other pathways. For instance, tillage variables are presented relative to no-till with manure surface-applied. More intensive tillage leads to more erosion, but also more effective incorporation of manure, i.e., higher P losses in soil and lower P losses from manure. The overall impact on P loss from reducing tillage is large, but because the P-index counts soil-bound P losses less, the impact on the P-index is smaller.

Comparison with basin-wide SWAT modeling

The results for our base scenario are similar, but not fully congruent with results from the SWAT modeling of the Lake Champlain basin conducted for the EPA's TMDL (US EPA, 2016). Differences with the Lake Champlain TMDL are difficult to fully compare since results are aggregated by lake segment (rather than sub-watershed) and most Vermont-side lake segments include areas outside the scope of our field-shapes dataset (which only covers four counties in Vermont). Additionally, the TMDL load allocations calculate P loading of the lake segments, while the P-index only calculates P delivered to waterways. Finally, base load calculations for the Lake Champlain TMDL were conducted for 2001–2010, while our model utilizes data on crop fields and practices for 2015.

However, comparison is possible for a few sub-watersheds or combinations of sub-watersheds highlighted in Table 5. Our source data estimate lower agricultural land use than the data used for the TMDL. Changing agricultural land use in the region only accounts for a small portion of this difference; the 2016 National Landcover Dataset shows only a small decrease in agricultural land cover compared to the 2011 version used in the TMDL analysis. On the other hand, agricultural land cover as reported in the Census of Agriculture declined by about 10% between 2012 and 2017 (NASS, 2021). Visual inspection of aerial imagery in comparison with our field delineations and the NASS land cover categorizations indicates that both datasets include some false positives and some false negatives for agricultural land cover.

Table 5. Comparison of study results with the SWAT model for the Lake Champlain TMDL

Only the sub-watersheds that are entirely within the study area are shown.

Our base model estimates for P loss per hectare of farmland are somewhat higher in aggregate to those in the TMDL, though there are significant discrepancies between watersheds. This shows that our base scenario, without NMP, is likely an over-estimate of total P losses. This makes sense, given that NMP was often practiced in Vermont even before regulatory changes made it mandatory for most farms. This also underlines our finding that targeted reductions may be difficult to achieve in some watersheds since our high reductions found from NMP are likely overestimates.

Discussion

Our results confirm the observation that P loading reduction targets for some sub-watersheds of the Lake Champlain basin may be difficult or impossible to achieve without a reduction in the agricultural footprint and/or size of the Vermont dairy herd. While maximum expected reductions in P loss are sufficient to meet reduction targets in some lake segments, they come up far short in others, especially Missisquoi Bay. Cattle numbers in the four county study area have been declining in recent years, falling by more than 14% between 2012 and 2017, accompanied by similar drops in cropland harvest for forage, and cropland harvest overall (NASS, 2021). However, a more intentional approach to herd size management may be warranted, particularly given a heterogeneous legacy of soil P accumulation.

Most previous spatial analyses of P loss on a watershed scale assume that soil phosphorus is distributed homogenously across the landscape. While this is an important simplifying assumption, our results suggest that the heterogeneity of legacy soil test P is an extremely influential factor in P loss from farm fields. Spatial analyses can easily detect ‘hotspots’ where soil type, rainfall, and proximity to waterways increase the risk of P loss. Additionally, we find that ‘invisible’ hotspots, caused by high legacy P in soils, may be just as significant for targeting conservation interventions.

Wironen, Bennett, and Erickson (2018) demonstrated large P surpluses throughout the Lake Champlain basin, driven in more recent years by imports of animal feed. Likewise, Ketterings, Kahabka, and Reid (2005) showed that nearly half of farm fields in New York have ‘high’ or ‘very-high’ soil phosphorus levels, with dairy-producing areas having the highest levels. High-soil P concentrations are generally caused by long-term repeated application of manure at high rates. Manure application history can be extremely variable between, and even within, fields due to variability in animal stocking rates and ease of access to fields for manure-spreading (Page et al., 2005). The soil test data used show ~90% of crop fields have soil test P levels at or below agronomic optimum levels, while a small fraction have values that are extremely excessive. Identification of sites of high legacy P may play a crucial role in mitigation efforts.

For example, Winchell et al. (2011) used a SWAT model to identify critical source areas for P loss in the Missisquoi Bay sub-watershed. Their work showed that targeting BMPs to the 20% most vulnerable fields could reduce P loss by 50–198% more than by applying these interventions randomly over the same number of fields. Our results suggest that even larger efficiency gains could be attained through such targeting. Their analysis with regard to ‘reduced manure P’ as a BMP may be an especially large underestimate. We find a 14% decrease in P loss from the ‘NMP’ scenario, which re-allocates manure P from fields with high P-index to lower P-index, based on the Vermont Required Agricultural Practices regulation.

Another issue deserving of more attention is the impact of pattern tile drainage commonly utilized on poorly and somewhat-poorly drained crop fields throughout the USA, including Vermont. The impacts of pattern tile drainage on nutrient losses are complex, not well-understood, and insufficiently handled by current modeling tools (Radcliffe et al., 2015; Wang et al., 2020). Pattern tile drainage reduces surface runoff and erosion, and thus the quantities of P lost through these pathways. However, tiles also create additional pathways for P loss through sub-surface drainage. Models generally assume that these losses are relatively small. For example, the Vermont P-index calculates subsurface drainage loss as 20% of the quantity of surface loss to the field edge. However, the literature is inconsistent with regard to the overall impacts on P loss (King et al., 2015).

In addition, important interactions between management and soil characteristics may complicate simplistic assessments. For example, most tile-drained crop fields in the Champlain basin are clay-textured and receive large applications of livestock manure. High clay fractions can lead to soil macropores (Beauchemin, Simard, and Cluis., 1998) which can serve as a direct path for manure-contaminated water into tile-drains, ditches, and streams (Dean and Foran, 1992; King et al., 2015; Shipitalo and Gibbs, 2000). While it was long assumed that most P loss occurred through overland flow, some measurement studies observe more than half of lost via sub-surface drainage (Smith et al., 2015).

Further complicating modeling efforts, some farming practices designed to attenuate P losses from erosion and surface runoff, such as conservation tillage and manure injection, may cause increases in P loss through subsurface drainage. Conservation tillage facilitates the formation of soil macropores, while manure injection places manure deeper into the soil profile. Within the study area, White et al. (2021) found that minimum tillage, cover crops, and manure injection on corn crops reduced P loss in surface runoff by 78%, which compares to a ~76% reduction in total P loss predicted by the coefficients in Table 4. On the other hand, P leaching into lysimeters increased by 121%. If all P leached below the root zone were conveyed into streams by pattern tile drainage, then this combination of BMPs would increase total P loss by 80%, rather than reducing it. Other watershed and plot-scale studies also show that conservation agricultural practices on drained cropland may decrease surface loss but increase sub-surface losses (Griffith et al., 2020; Jarvie et al., 2017; Smith et al., 2015). This problem is also reflected in the TMDL BMP scenario analysis tool developed using SWAT (Tetra Tech Inc., 2015), which predicts a 48% reduction in P losses from the BMPs used in that study.

In our study, the Vermont P-index predicts slightly lower (~−1%) total P-index values and slightly higher total P lost (~16%) from crop fields with tile drainage in the base case, all else held equal. This increase reflects the balance of reduced erosion and runoff from improved drainage, but that rainfall near tile drain lines can directly convey particles into the drain lines. The Vermont P-index currently estimates that this occurs over only 20% of the field. The current version of the Vermont P-index is not capable of handling interactions between conservation agriculture practices and pattern tile drainage.

Conclusion

In this study, we scaled-up a field-level model to examine farm P losses on a watershed scale. Results show that most sub-watersheds in the Vermont side of the Lake Champlain basin are unlikely to achieve their targeted reductions without reducing livestock numbers and the extent of agricultural land, or utilizing novel, yet-developed BMPs. We further show that spatial heterogeneity in soil legacy phosphorus may be an underexamined lever for targeting interventions, and especially for targeting land retirement or shifting to low-intensity cropping systems.

Simulating site-assessment tools on a large-spatial scale is, to our knowledge, a novel approach to examining nutrient loss from farmland at the basin scale. Because each regional P-index is slightly different, our simulation code cannot be directly applied in other places, but our general approach might be. Similarly, it could be utilized for site assessment tools including models of nitrogen loss and carbon sequestration on farmland.

Future research using this approach could be helpful for validating and/or improving site assessment tools. Furthermore, this work could be cross-validated with NMP records which utilize these same site assessment tools. Our framework could be helpful for estimating total reductions and financial costs of performance-based programs for reducing nutrient losses and for targeting interventions to highest-impact locations. Lastly, it provides a flexible framework for examining how large the impacts of proposed interventions would need to be to make these interventions meaningful.