Investigating tradeoffs between competing agricultural objectives using multi-objective optimization under uncertainty

Management of groundwater quality in agricultural areas requires tradeoffs between competing objectives. These objectives include economic benefit, respect for regulatory-imposed water quality standards, and an equitable distribution of revenue that is foregone in order to meet these standards. It also requires recognition of the fact that the calculated outcomes of a proposed management strategy are accompanied by uncertainty. In this study, we apply formal multi-objective optimization under uncertainty to investigate these tradeoffs as they apply to nitrate-N leaching rates in the northern Pineland Sand Aquifer System in Central Minnesota. This region is important both to agricultural producers and the Anishinaabe Peoples of White Earth Nation. The methodologies that we demonstrate exemplify a scientific framework in which simulation-based, management-pertinent analyses can support negotiations between different stakeholder groups wherein the risks and consequences associated with different management strategies are clearly delineated. At the same time, it is able to advance stakeholder understanding of the environmental factors that are salient to optimal groundwater management, and the necessity to accept uncertainty, and hence risk, when exploring management options. It also provides insights into sources of uncertainty, and strategies that can most effectively reduce it.

1 Introduction

Nitrate-Nitrogen (NO₃-N) contamination of groundwater systems is a major concern in agriculturally dominated watersheds. In the rural Midwest of the United States, fertilizer applications commonly result in NO₃-N, carried by infiltrating water, leaching into groundwater systems (Cartwright et al., 1991; Novotny, 1999). The US Environmental Protection Agency (US EPA) sets the maximum NO₃-N concentration level for drinking water at 10 mg L^-1 in order to reduce impact on human health from conditions such as methemoglobinemia, also known as blue baby syndrome (United States Environmental Protection Agency, 2007). Ward et al. (2018) conducted a literature review on current research into the human health impacts of NO₃-N exposure. In addition to methemoglobinemia, this study documented identification of possible linkages between groundwater NO₃-N and (among other things) an increased risk of colon and colorectal cancer (Schullehner et al., 2018; McElroy et al., 2008; Espejo-Herrera et al., 2016; De Roos et al., 2003), thyroid disease (Ward et al., 2010; Tajátková et al., 2006; Aschebrook-Kilfoy et al., 2012; Gatseva and Argirova, 2008), and neural tube defects (Brender et al., 2004; Brender et al., 2013; Croen et al., 2001; Dorsch et al., 1984; Arbuckle et al., 1988).

NO₃-N in groundwater can discharge into lakes, rivers, and streams. NO₃-N concentrations in rivers and streams are often elevated during dry seasons when their flows are baseflow dominated. This can have adverse ecological consequences. Kincheloe et al. (1979) found that NO₃-N is mildly toxic after 30 days at concentrations of 1.1–4.5 mg NO₃-N L^-1 for rainbow trout, chinook salmon, Lahontan cutthroat trout, and steelhead trout at the developing egg to early fry stages. Their findings listed a NO₃-N concentration of 2 mg L^-1 as limiting reproductive success of some salmonid fish. Chronic exposure to NO₃-N concentrations in excess of 6.25 mg L^-1 have been observed to impact fry weight of lake trout (McGurk et al., 2006). The State of Minnesota has set a draft ecological chronic exposure standard of 5 mg NO₃-N L^-1 for salmonid waters and 8 mg NO₃-N L^-1 for non-salmonid waters (Monson, 2022). At the time of writing, a standard has yet to be set by the Minnesota Pollution Control Agency (MPCA).

The study that is documented herein examines the impact of land management on NO₃-N concentrations in the Pineland Sands aquifer system in the north-central region of Minnesota. See Figure 1. This is a glacial system consisting of sand and gravel aquifers and adjacent glacial diamicton aquitards (Bauer et al., 2016; Lusardi, 2016; Lusardi, 2018). The aquifer system feeds the headwaters of the Upper Mississippi River. Increasing NO₃-N concentrations in ground and surface waters pose risks to human and ecological health. They also have adverse cultural repercussions for the Anishinaabe Peoples of the White Earth Nation who enjoy federally protected hunting, fishing, and gathering rights within the area, and who use the lakes, rivers, and streams for harvesting of wild rice, a crop that is intrinsically linked to Anishinaabe spiritual practices and culture (Katanski, 2017).

Figure 1

Figure 1. Map of study area with respect to the State of Minnesota.

The Pinelands Sands aquifer system is overlain by sandy soils that are favorable to farming. Since the 1970s, the area has experienced a continuous conversion of land use from dry-land and forest to irrigated agriculture. Recent years have seen an acceleration in this conversion. For example, the Straight River Watershed (see Figure 1) has experienced an approximately 84% increase in irrigated acreage from 1992 to 2021. Meanwhile, other parts of the study area are being considered for future agricultural land use. Common row crops in the region include potato, corn, soybean, dry beans, and spring wheat, with other crops being grown to a lesser extent (Stark et al., 1994; Stroom, 2024).

NO₃-N concentrations are regularly monitored in streams and in some groundwater wells. Monitoring in the Straight River Watershed has revealed increasing NO₃-N concentrations in both ground and surface waters. In 2009, measured NO₃-N concentrations exceeded the human health standard of 10 mg L^-1 in a well that supplies the City of Park Rapids. As a result, a deeper well was drilled (Stroom, 2024). Because of rising public concerns over increased NO₃-N concentrations both within the study area and more broadly, the State of Minnesota has enacted a nutrient reduction strategy. This strategy sets a goal of a 45% reduction of nitrogen load into the Mississippi River by 2040 (MPCA, 2014).

In common with environmental stewardship in general, management of NO₃-N levels in ground and surface waters presents a multi-faceted technical and social dilemma. On the one hand, economic benefits flow from increased farming, and from increased use of fertilizers to support growth of high demand crops such as potatoes. On the other hand, this cannot be allowed to compromise drinking water quality, nor the ecological health of high value waterways and the ecosystems that they support, particularly where these ecosystems have high cultural significance. Management of a system such as the Pinelands Sands aquifer therefore implies a tension between competing objectives. Management is further complicated by the fact that different societal groups may place different weights on some of these objectives.

Negotiations and policy development that accompany groundwater management can be facilitated if certain key elements are established as the premise for these discussions. The first is an acceptance that management of the system be science-based. This requires mathematical processing of data to enable harvesting of management-pertinent information from these data. Information that is acquired in this way is expressed in predictions of system behavior under different management options in order that the consequences of these options can be assessed by all interested parties. Just as importantly, science-based data processing should delineate the management repercussions of information insufficiency. This is generally expressed as predictive uncertainty. Whether uncertainty is associated with particular management outcomes, or with comparative management outcomes, its quantification is essential to assessment of management risks. We use the word “risk” to loosely denote the probability of occurrence of an unwanted management outcome multiplied by the cost of that outcome. Obviously, the association of cost with an outcome is value-dependent; it may differ for different management stakeholders.

A second requirement for facilitation of discourse between management stakeholders is development and acceptance of an intellectual framework for these discussions. Ideally, at the core of this framework are scientifically-based predictions of management outcomes. Also central to this framework is the notion of tradeoff, and the ability to assess the implications of tradeoffs by different stakeholder groups. A tradeoff may be required between pursuit of one management outcome and constraints that may be emplaced on that outcome. A tradeoff may also be required between tolerance of risk and pursuit of a certain management goal. Negotiations between interested parties require that the costs and benefits of any management goal be defined, that tradeoffs be analyzed, and that the risks associated with tradeoffs be clearly enunciated.

It is acknowledged that decisions are rarely made by selecting a point on a tradeoff curve. Instead, they are based on a “narrative” (Kay and King, 2020). Nevertheless, the development of a narrative that embodies scientifically-based collection and processing of data while acknowledging limitations in information that can be extracted from these data, can provide a meeting ground where those with different values can negotiate in good faith, with clear separation of those values from prognostications of future system behavior.

In this paper, we present a methodology for assessing the current state of a groundwater system with regard to NO₃-N contamination, and of quantifying uncertainties associated with that assessment. We then address the issue of how management of that system can attempt to reconcile the interests of different parties who are stakeholders in that system, by displaying the benefits, costs, and tradeoffs that its management requires. These methods are demonstrated using a groundwater model of the northern Pinelands Sands aquifer. However, the methodology is general and is easily applied elsewhere. The intention of this paper is to demonstrate an intellectual framework which can help further provide a better distillation of the scientific elements that are integral to the NO₃-N management decision-making process. This process typically requires significant technical and social input and is challenged by the social ramifications from affected growers. Therefore, this study focuses on the technical input which this framework can help provide, rather than implementing a decision-making process.

The paper is organized as follows. First, we briefly describe technologies that are deployed in this study. Two of these are particularly salient to environmental (particularly groundwater) management support. The first of these is an iterative ensemble smoother. This provides an efficient mechanism for extraction of information from field measurements of groundwater behavior through simulator history-matching. The second is multi-objective optimization under uncertainty. This provides a means of examining tradeoffs between competing management objectives, and between pursuit of those objectives and acceptance of risk. We then use these methods to examine costs, benefits, and constraints that are associated with management of the northern Pinelands Sands Aquifer. This is followed by a discussion of the results of that analysis as it pertains to the Pineland Sands Aquifer itself, and aquifer management in general. Conclusions follow that.

None of the simulation or data-processing methods that we discuss in this paper are new. Hence descriptions of them are brief. However, as far as the authors are aware, the manner and location of their deployment, and the issues that they address, have not been published in the scientific literature. Furthermore, in this paper we take the opportunity to present them to a different audience from that in which they are commonly discussed in order (we hope) to encourage their wider and more fruitful adoption. In doing this, we attempt to emphasize their utility in providing an intellectual framework that enables discussion amongst parties whose interests and values may differ, but who inhabit the same area while sharing the resources that this area offers.

2 Methodologies

2.1 General

Groundwater simulation is used ubiquitously to provide scientifically-based support for environmental management. However, its use is often contentious. Different sides of an argument often deploy different groundwater models that make predictions which are favorable to them, based on the same data. Faced with the spectacle of dueling models, a decision stakeholder cannot be blamed for treating groundwater modelling with skepticism.

Doherty and Simmons (2013) remark that this situation can be rectified through adherence to appropriate decision-support groundwater modelling metrics. These are based on the premise that predictions of future groundwater behavior cannot be made with certainty. The task of decision-support groundwater modelling is therefore to quantify the uncertainties of management-salient groundwater predictions while reducing them as much as available data allows. Groundwater modelling, as a decision-support activity, fails if its uncertainty intervals are so narrow as to preclude outcomes that are indeed possible. On the other hand, groundwater modelling is of limited decision-support utility if uncertainty estimates are artificially broad.

Uncertainty is reduced through acquisition of information. The task of groundwater modelling is therefore to harvest information, and to deliver this information to decision-makers in a form that they can use. At the same time, decision-makers must be made aware of the consequences of information insufficiency. The latter is expressed by the uncertainties that remain in predictions of management interest after all available data has been processed.

Within a groundwater model, information is carried by its parameters. Parameters are, by definition, adjustable. Their prior probability distribution expresses what is known, and what is unknown, about system hydraulic properties and boundary conditions as an outcome of conceptual model development and site characterization. More information is then harvested through history-matching. Model parameters are thereby constrained by the necessity for pertinent model outputs to replicate measurements of system behavior (to the level of noise that is associated with those measurements). Their post-history-matching probability distribution is referred to as their posterior probability distribution.

2.2 Simulation

In the present study we employ MODFLOW-NWT (Niswonger et al., 2011) to simulate movement of water and MT3D-USGS (Bedekar et al., 2016) to simulate transport of NO₃-N. The use of MODFLOW/MT3DMS, and programs that have descended from it, in particular, MODFLOW-USG (Panday et al., 2013) and MODFLOW 6 (Langevin et al., 2017) to simulate movement of NO₃-N is relatively commonplace in the groundwater literature. See, for example, Daniel et al. (2023), Yang et al. (2021), Llopis-Albert et al. (2014), Reed et al. (2018), Shakeri et al. (2023), Sieczka et al. (2018), Almasri and Kaluarachchi (2007), Shamrukh et al. (2001), Ameur et al. (2021), and Correa-González et al. (2023).

2.3 History-matching

The term “history-matching” describes the adjustment of model parameters in order to achieve a good fit between field measurements and their model-calculated counterparts. This is often loosely referred to as “model calibration”. However, “calibration” refers to only one type of history-matching. This is the search for a unique parameter field for which (theoretically at least) all model predictions are of minimum error variance. That is to say, these predictions have minimized potential for predictive error. This is because they lie at the centers of their respective posterior probability distributions. Possible predictive errors (which may be considerable) are therefore symmetric with respect to these predictions. See Doherty (2015) for more details.

In the present study we undertake Bayesian history-matching by attempting to sample the posterior probability distribution of model parameters. When a prediction is made using a model that is history-matched in this way, it is then made using all of these samples. The posterior predictive uncertainty of that prediction is thereby sampled. Bayesian history-matching is implemented using the PESTPP-IES iterative ensemble smoother; see PEST Development Team (2025), White (2018) and Chen and Oliver (2013) for algorithmic details. To initiate the process, a modeler provides samples of the prior parameter probability distribution. A few hundred samples (i.e., realizations) are normally enough; however, where data are plentiful, more samples than this may be required to fit the history-matching dataset and to quantify parameter uncertainty. Through iterative adjustment, these samples of the prior parameter probability distribution are morphed into samples of the posterior parameter probability distribution. All of these samples allow the model to adequately replicate field measurements.

Two major advantages accompany use of an iterative ensemble smoother (IES) for groundwater model history-matching. The first is its ability to manipulate many parameters (up to tens of thousands). The second is the small numerical burden that is required for this manipulation. The number of model runs required per iteration of the history-matching process is equal to the number of parameter realizations that require manipulation. As stated above, this is normally of the order of a few hundred. The number of iterations required for attainment of an acceptable fit with a history-matching dataset varies from case to case. In the present study, three iterations were sufficient.

The endowment of a groundwater model with many parameters is fundamental to the achievement of its decision-support potential. There are two reasons for this. The first is that it relieves a modeler of the necessity of making unrealistic assumptions pertaining to patterns of subsurface heterogeneity such as that of piecewise homogeneity that is often employed as a basis for parsimonious parameterization. Instead, heterogeneity can “find itself” during the history-matching process subject to regularization controls that require realistic expression of this heterogeneity. The second arises from the nature of the history-matching task itself. When quantifying parameter and/or predictive uncertainty, representation in a model of the heterogeneity that cannot be uniquely estimated is just as important as representation of that which can. The former is required to express history-matching nonuniqueness; this is normally the major contributor to groundwater model parameter and predictive uncertainty. As stated above, the numerical cost of accommodating a large number of parameters is small when using IES.

Conceptually, a Bayesian approach to decision-support groundwater modelling could invite more stakeholder participation than deterministic approaches used in the past that were based on predictions made using a calibrated model (in which uniqueness was often pursued through parameter parsimony) whose uncertainties were considered to be either small, or irrelevant to decision-making. This could be achieved, for example, through collegiate setting of prior parameter probability distributions, and/or in testing the compatibility of an unwanted hypothesized management outcome with that of observed past system behavior. This, it is hoped, could eliminate the inglorious specter of competing models.

2.4 Constrained optimization

The task of optimization is to adjust so-called “decision variables” in order to maximize or minimize a quantity of interest. In optimization parlance, calculation of this quantity requires evaluation of an objective function. Where the optimization process is constrained, then the values assigned to decision variables must simultaneously ensure that inequality relationships for user-specified dependent variables are maintained. In environmental management contexts, costs may be minimized, or benefits may be maximized. Meanwhile, constraints pertain to the values of system states or fluxes that, if exceeded, can be detrimental to environmental and/or human health. In the example below, we maximize farm revenues while simultaneously ensuring that NO₃-N levels at certain observation points are lower than regulatory thresholds.

A quantity that is undergoing maximization/minimization may not require a complex model for its evaluation. It may, for example, be a simple summation of items that contribute to benefit or cost. In contrast, respect for constraints often requires the use of an environmental model. The model may calculate, for example, the value of an environmental variable at one or a number of measurement sites; in the case presented below this is NO₃-N concentration. The constrained optimization process must then guarantee that these values are not exceeded as the objective function is optimized.

So-called “global” optimization methods are often deployed where relationships between dependent and decision variables are highly nonlinear, and/or where the existence of local objective function maxima/minima complicate the optimization process. Global methods gain their robustness through strategic use of randomization to probe decision variable space. Different global methods differ in their implementation of randomization and in how they ensure that learning is maximized as the optimization process progresses. Unfortunately, their model run cost is considerably greater than that of commonly used methods such as sequential linear programming (SLP) and sequential quadratic programming (SQP) that perform well in less demanding optimization contexts.

An advantage of global optimization methods is that they can be readily adapted to the exploration of tradeoffs between pursuit of different objectives, and/or optimization of an objective function and respect for constraints. They do this by defining a Pareto front. The dimensionality of this front increases with the number of items for which a tradeoff is sought.

Suppose that a global optimization process is exploring the tradeoff between maximization of two objective functions. Under these conditions, the Pareto front is the locus of points in decision space for which it is not possible to raise the value of one objective function without simultaneously lowering the value of the other. Ideally, the “best” set of decision variables lies somewhere along this front, with “best” reflecting relativity of values ascribed to the two objective functions. Where values ascribed to these objectives by different stakeholder groups are different, then exploration of the Pareto front allows decision participants to view for themselves the extent to which attainment of one objective impedes attainment of a competing objective.

In the study that is discussed herein we employ PESTPP-MOU (PEST Development Team, 2025) to perform global optimization; see White et al. (2022) for algorithmic details. This provides a number of global optimization methods. The Particle Swarm Optimization generator option (Eberhart and Kenedy, 1995) is used in the present study, together with the Non-Dominated Sorting Genetic Algorithm II (NSGA-II) (Deb et al., 2002) for exploration of the Pareto front. 74 iterations were used and a population or “particle” size of 750 was chosen for this experiment. Each particle represents a set of decision variables. Decision variables associated with each particle are iteratively improved as the optimization process progresses. The NSGA-II algorithm spreads these particles along an evolving Pareto front as evenly as possible. The greater the number of “particles”, the greater are the run time requirements of the optimization process. However, the greater is the surety that the Pareto front is properly characterized.

2.5 Constrained optimization under uncertainty

Respect for uncertainty in optimizing management of groundwater systems has been a subject of active research interest for many years. See, for example, Gorelick (1997), Bayer et al. (2010), Coulon et al. (2022), Coulon et al. (2024), White et al. (2020), Fienen et al. (2024) and papers cited therein.

Where a model’s parameters are uncertain, quantities that are calculated by the model inherit this uncertainty. These include quantities that are used to constrain an optimization process. Enforcement of optimization constraints then becomes probabilistic. That is to say, a user must specify his/her tolerance for a constraint being violated as an objective function is maximized or minimized. The constraint is then respected to this level of tolerance. To simplify the following discussion, we assume that constraints are respected if model-calculated quantities are less than complementary user-specified thresholds.

To see how this works, suppose that a model is endowed with N_p parameter fields, these being the outcomes of Bayesian history matching. Hence all of these N_p parameter fields are samples of the posterior parameter probability distribution. As such, they are “realistic” (according to the prior parameter probability distribution) whilst also being “behavioral” (that is, they support a satisfactory fit between pertinent model outputs and measurements of historical system behavior). Suppose also that a decision-maker adopts a 100% risk-averse attitude to constraint enforcement. Then, when the optimizer runs the model based on a certain choice of decision variables in order to evaluate model-calculated values of constraints, it runs the model N_p times; on each occasion that it runs the model, it employs a different posterior parameter field. The model-evaluated value of the constraint must be less than the constraint threshold for all N_p parameter fields for the constraint to be decreed as respected. The same applies to all other constraints.

Suppose now that a decision-maker adopts 95% risk averse attitude to constraint enforcement. Then the constraint is decreed to be respected if pertinent model-evaluated quantities are less than constraint thresholds for 95% of the N_p model runs. The same applies to other levels of risk aversion. Risk neutrality implies a 50% risk aversion stance.

The concept of tradeoff, and its exploration through definition of a Pareto front, can be applied to risk. This is because a high level of risk aversion generally inhibits the extent to which an objective function can be maximized or minimized, while risk neutrality or risk acceptance provides more opportunities for decision variable adjustment in order to optimize an objective function. In the study that is documented below, we decree the risk stance of a decision maker to vary between 0.0 and 1.0, with 1.0 corresponding to 100% risk averse.

Because the running of a model during a constrained optimization-under-uncertainty process becomes the running of N_p models, where N_p is the number of parameter fields that sample the posterior parameter probability distribution, and because the model must be run many times using many different values of decision variables in order to maximize/minimize an objective function and/or explore a Pareto front, optimization under uncertainty is an extremely numerically intensive process, especially when implemented using a global optimization method. PESTPP-MOU provides a number of options for reducing the numerical burden of this process. These include re-use of uncertainty interval widths in successive iterations of the optimization process, with intermittent updating of these widths at user-specified iteration intervals. See PESTPP-MOU documentation for further details.

3 Model construction and history-matching

3.1 General

Concerns over increasing NO₃-N contamination of the Pinelands Sands aquifer system are discussed in the introduction to this paper. The focus of modelling work that is described herein is the northern Pineland Sands region. Figure 1 shows the boundary of the model domain. Within the study area is Pine Point township which is at the southeasternmost corner of the White Earth Nation Reservation; the White Earth Nation are partners in this investigation. Because of their cultural linkages to lakes and streams within the study area, they are also groundwater management stakeholders. So too are farming communities whose activities are impacting groundwater while providing them with economic benefits, as are state and local authorities who are tasked with managing land and water resources for the region. Modelling that is discussed herein is targeted at exploration of management options, and the tradeoffs that implementation of these options require.

3.2 Simulation specifications

Figures 5,6 provide an overview of the methodology used in this study. As mentioned above, groundwater flow is simulated using MODFLOW-NWT while solute transport is simulated using MT3D-USGS. The grid for both of these simulators is structured, with a uniform cell size of 200 m × 200 m. It is comprised of three layers. The three layers represent composites of upper, middle, and lower glacial formations based on approximate age of the glacial formations that collectively comprise the Pineland Sands Aquifer System. A similar approach to model construction and layering was adopted by Fienen et al. (2022) in a USGS study of the Wisconsin Central Sands area. The total number of active cells is 115,102. See Supplementary Material for more information.

Groundwater flow is considered to be steady state. This assumption is unlikely to adversely affect model-evaluated linkages between regional sources of NO₃-N and the areal disposition of NO₃-N concentrations within the aquifer as regional flow directions do not change greatly from season to season. Groundwater recharge used in the model comes from the Minnesota Soil Water Balance Model (Smith and Westenbroek, 2015) run over 2002–2015. The River Package (RIV) is used to represent lakes, rivers, and streams in the model area. The General Head Boundary Package (GHB) is used to represent outer model boundary conditions. Hydraulic heads for GHB boundary conditions are gleaned from the County Geologic Atlas groundwater contours for Hubbard, Becker, and Wadena counties (Bradt, 2023; Budde, 2024; Budde and Baratta-Person, 2024); in general, these heads pertain to distances of 1 km from outer model boundary cells. A small number of cells in layer 1 (36 cells) employ the Drain Package (DRN) for outer boundary conditions where contoured hydraulic heads are close to cell bottom elevations.). The multi-node well package (MNW2) is used to represent high-capacity pumping wells. The governing equation for 3-D steady state groundwater flow is defined in Equation 1.

$\frac{\partial }{\partial x}\left(K_{xx}\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yy}\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_{zz}\frac{\partial h}{\partial z}\right)+W=0(1)$

In this equation, K_xx, K_yy, and K_zz are the hydraulic conductivities in the x, y, and z directions, h is the hydraulic head, and W is the volumetric flux per unit volume from sources/sinks (Harbaugh, 2005).

MT3D-USGS simulates NO₃-N transport over a period of 100 years. This is long enough for NO₃-N concentrations within the model domain to stabilize and hence achieve steady state conditions. Equation 2 describes the governing equation for 3-D advective-dispersive-reactive solute transport in groundwater flow systems.

$\frac{\partial }{\partial x_i}\left({\theta D}_{ij}\frac{\partial C^k}{\partial x_j}\right)-\frac{\partial }{\partial x_i}\left(\theta v_i{C}^k\right)+q_s{C}_s^k+\sum R_n=\frac{d\left(\theta C^k\right)}{\partial t}(2)$

where $\theta$ is porosity, $D_{ij}$ is the hydrodynamic dispersion coefficient tensor, $C^k$ is species k’s dissolved concentration, $x_i$ and $x_j$ are the distances along the respective coordinate axis, $v_i$ is the average linear velocity, $q_s$ is the volumetric flux per unit volume of sources/sinks, $C_s^k$ is the concentration of the source/sink flux for species k, $\sum R_n$ is the chemical reaction term, and t is time (Zheng and Wang, 1999).

Advection is simulated using the standard finite difference method. This option is faster than alternative solution devices offered by MT3D-USGS, an important consideration given the centrality of optimization to the present study and the many model runs that this requires. Denitrification in this study was treated as an uncertain, spatially variable parameter represented by 1,559 pilot points for each layer spaced 1 km apart. The denitrification rate for each pilot point is centered on a first-order rate constant of 7.26495 × 10⁻⁵ d^-1, representing the median denitrification rate calculated from field data measured by Puckett and Cowdery (2002) in the surficial outwash aquifer within the Otter Tail Watershed southwest of the model area. This was assumed to be a representative denitrification rate due to the prevalence of outwash as the surficial sediment for much of the model domain including the areas of cultivated cropland which are the focus of this study.

The value assigned to each pilot point was allowed to vary between a lower bound of 5.44871 × 10⁻⁵ d^-1 and an upper bound of 9.08119 × 10⁻⁵ d^-1. This range reflects an uncertainty interval of ±25% from the median value measured in the field, which was chosen to represent upscaled regional heterogeneity for the composite layers that form the model domain. Local patches of higher and lower denitrification were observed in Puckett and Cowdery (2002). However, given the high degree of variability both vertically and horizontally over short distances within the same aquifer, these rates were assumed to not be representative of the scale with which is being modeled in this study.

Longitudinal dispersivity is set to 10% of cell length; horizontal transverse dispersivity is set to 10% of this value and vertical transverse dispersivity is set to 1% of this value. The effective molecular diffusion coefficient is assigned a value of 5 × 10⁻⁵ m² d^-1 (Frind et al., 1990). Porosity is set to 0.25 throughout the model domain. It can be argued that (except for the molecular diffusion coefficient) these model inputs should be declared as adjustable as they can affect model outcomes (because of denitrification). However, this is not done in the present case as spatial variability of these quantities is considered to have a secondary influence on the areal disposition of model-evaluated NO₃-N concentration. Likewise, dispersivity in numerical modeling of advective transport typically mimics the action of unrepresented heterogeneity in aquifer hydraulic properties (Zheng and Bennett, 2002). Because lateral variations of hydraulic conductivity heterogeneity are specifically represented in the model, an assumption of uniform dispersivity avoided overrepresentation of the dispersive component of contaminant transport. Nonetheless, for reasons discussed above, methodologies that are discussed herein would not suffer an undue numerical burden if these model inputs’ uncertainties were considered.

NO₃-N is introduced to 7,167 model cells containing cultivated cropland using the SSM package of MT3D-USGS. This land use type is identified using the NLCD 2019 dataset for land cover in Minnesota (Dewitz, 2021); see Figure 2. In the present study, NO₃-N leaching rate is a decision variable, and is therefore adjusted during constrained optimization. Parameterization of NO₃-N leaching rate is discussed below.

Figure 2

Figure 2. Map of cultivated cropland taken from NLCD 2019 dataset. Locations of the NO₃-N leaching rate decision variable pilot points as well as the groundwater quality constraints are also shown in this figure.

3.3 History-matching

The history-matching dataset for the northern Pineland Sands model is comprised of time-averaged heads in 541 wells, together with two measurements of stream baseflow. Hydraulic head observations were obtained from the Minnesota Well Index and Minnesota Cooperative Groundwater Monitoring Network (MDH, 2025; Minnesota Department of Natural Resources, 2025). Baseflow is obtained by applying a Lyne-Hollick baseflow separation filter to streamflow hydrographs taken from the Minnesota Cooperative Stream Gaging Network Database (Lyne and Hollick, 1979; Minnesota Department of Natural Resources and U.S. Geological Survey, 2025). The latter are matched to MODFLOW-calculated inflows into RIV package boundary conditions in cells upstream from pertinent stream gauging stations. Observation well and stream gauge locations are depicted in Figure 3. In this figure, wells are colored according to the model layer in which the majority of their well screen resides.

Figure 3

Figure 3. Locations of observation wells and stream gaging stations that comprise the calibration dataset.

Model parameters that are adjusted in order for model outputs to match these measurements include horizontal hydraulic conductivities and vertical anisotropies in each of the three model layers. In all cases, pertinent hydraulic property values are ascribed to a suite of pilot points that are distributed throughout the model domain. Each pilot point hosts a horizontal hydraulic conductivity and vertical anisotropy parameter. 243, 214, and 181 pilot points are assigned to layers 1, 2, and 3 respectively. Prior pilot point horizontal hydraulic conductivities are obtained through a weighted average of sand and till hydraulic conductivities based on the fraction of sand out of total material at each points location. Their locations allow them to express spatial variability of sand content. See Supplementary Material for more information. Hydraulic properties assigned to pilot points undergo spatial interpolation to individual model cells. A total of 1,276 model parameters are ascribed to pilot points.

As is discussed above, a large number of model parameters allow representation of patterns of spatial heterogeneity. On the one hand, this large number supports attainment of a high level of fit between field measurements and pertinent model outcomes. On the other hand, it far exceeds the number of parameters that can be estimated uniquely. These pilot point parameters can therefore support stochastic representation of hydraulic property details to which predictions of management interest may be sensitive. These predictions are NO₃-N concentrations that are employed in the constrained optimization process that is described below.

Note that groundwater NO₃-N concentrations in monitoring wells, though available, were not included in the history-matching dataset. Conceptually, these are high information content, particularly with respect to denitrification rates. The reason for this is that the historical nitrate “source term” requires knowledge of past farming areas and practices with their accompanying leaching rates. This information is unavailable.

Another 26 model parameters are devoted to representation of GHB conductance, while one parameter is devoted to representation of DRN conductance. Further information on these parameters can be found in the Supplementary Material for this paper.

Definition of prior parameter probability distributions requires that all parameters be provided with mean values and standard deviations. These are listed in Table 1. Being spatially distributed, pilot points also require a covariance matrix. For all families of pilot points, this is based on an exponential decay of correlation with distance; the decay constant is 1.5 times the local average pilot point separation.

Table 1

Table 1. Prior means and standard deviations of log(to base 10) parameter values.

As discussed in Section 2, history-matching is implemented using the PESTPP-IES ensemble smoother. 500 realizations of the prior parameter probability distribution are adjusted in order to sample the posterior parameter probability distribution, resulting in 475 realizations of the posterior (Attrition of realizations is not uncommon when using ensemble methods; some random parameter fields challenge iterative solution procedures used by numerical models). Figure 4 shows four realizations of posterior horizontal hydraulic conductivity in layer 1.

Figure 4

Figure 4. Some horizontal hydraulic conductivity fields emerging from PEST++ IES history matching process. Four realizations of layer 1 horizontal hydraulic conductivity fields are shown: (A) base realization (an approximation to mean over all parameter fields), (B) realization 0, (C) realization 238, (D) realization 475. Many of the hydraulic property details are inherited from local variation in sand content.

Figure 5

Figure 5. Flow chart depicting the steps in the model building and uncertainty analysis process.

Figure 6

Figure 6. Flow chart depicting the steps in the optimization and multi-criteria analysis process.

Simulated water table elevations in the base realization obtained from history-matching show a dominant groundwater flow direction towards the Crow Wing River, which is the regional discharge outlet for the aquifer system. Local groundwater discharge is observed in the smaller rivers and streams within the model area such as the Straight River and Shell River. For most of the model domain which resides to the west of the Crow Wing River, groundwater flow directions are generally from the northwest towards the southeast of the model domain. See Supplementary Material for realizations of hydraulic properties in other model layers, model-to-measurement misfit for the groundwater heads component of the history-matching dataset, and the simulated water table elevations.

4 Optimization under uncertainty

4.1 Economic objectives and decision variables

We now formulate an optimization problem that represents the competing objectives that characterize management of the study area (and many other areas). On the one hand, there is an economic imperative to maximize farm productivity, and with it the prosperity of the area as a whole. On the other hand, because higher rates of farm productivity require greater application of nitrogen fertilizer, farm productivity is constrained by the necessity to preserve groundwater quality and the health of ecosystems that are fed by groundwater.

In order to simplify formulation of the optimization problem, we make the following approximations. First, we use NO₃-N leaching rates as a surrogate for farm revenue. By increasing the NO₃-N leaching rate in any part of the model domain, farm revenue in that area is thereby construed to be increased. In the constrained optimization process, NO₃-N leaching rates are formulated as continuous, but spatially variant, decision variables. Leaching rates are assigned to an array of 271 pilot points dispersed over pertinent parts of the model domain; see Figure 2. The leaching rate associated with each point is spatially interpolated to model cells in which crops are grown. NO₃-N is then introduced to the groundwater system using the SSM package of MT3D-USGS. It is important to note, however, that the leaching rate assigned to any pilot point is limited to a maximum of 150 kg ha^-1 yr^-1 during the optimization process. This corresponds to the NO₃-N leaching rate measured by Errebhi et al. (1998) for potatoes in Becker County, Minnesota based on the temporal fertilizer application strategy which produces the largest potato yield. The link between other leaching rates and crop types is discussed shortly.

The objective function to be maximized is defined in Equations 3a,b as follows:

${\mathrm{\Phi }}_N={\displaystyle \sum _{i=1}^{M_d}}{N}_i(\mathrm{3a})$

with the range of each N_i restricted such that:

$N_i\le 150;1\le i\le M_d(\mathrm{3b})$

In these equations M_d is the number of pilot points with which decision variables are associated (i.e. 271) while N_i is the leaching rate associated with the ith pilot point.

4.2 Constraints

Balanced against the objective of maximizing NO₃-N leaching rate (as a surrogate for farm revenue) are constraints on groundwater quality. NO₃-N concentrations are calculated by the model. For one set of leaching rates, these calculations are repeated 475 times using the 475 samples of the posterior parameter probability distribution obtained during PESTPP-IES history-matching to obtain uncertainty interval widths for the constraints. These widths are recalculated every five iterations. Constraints on groundwater NO₃-N concentrations are set at the 34 wells that are depicted in Figure 2. A further constraint is applied to the concentration of NO₃-N in groundwater that emerges into the Straight River reach that is also depicted in Figure 2. The Straight River was chosen as the only waterbody constraint imposed in the study area because of its importance for the region as a result of its salmonid water classification. Constraints were not applied to lakes in this demonstration study as this would require better simulation of factors which govern denitrification in these surface water bodies. For similar reasons, other streams and rivers in the region were not selected as constraints due to their non-salmonid water classification and/or because of their connection to flow-through lakes. In the present study, imposition of Straight River constraint assumes NO₃-N concentrations in lake outflow are zero due to denitrification in the lake. The constraint is calculated by dividing the modeled total mass flux for the entire lower reach seen in Figure 2 by total groundwater reach leakage into the Straight River and Straight Lake.

For flow into the river reach, the constraint is set to 5 mg L^-1, as this reach is classified as salmonid water. In the optimization processes that are reported below, two sets of constraints are applied to NO₃-N concentrations at the 34 groundwater wells, namely, 5.4 mg L^-1 and 8.0 mg L^-1. The former is denoted as Mitigation Level 1 by the Minnesota Department of Agriculture (2019) while the latter is denoted as Mitigation Level 2. These are levels at which voluntary NO₃-N mitigation procedures are requested. If Level 2 practices meet with no success, adoption of appropriate mitigation practices becomes a regulatory requirement.

Constraints can be formulated using Equation 4.

$\left[N_j\right]\le \left[N_j^T\right];j\le 1\le M_m(4)$

where [N_j] is the concentration of NO₃-N at the jth monitoring point, [N^T_j] is the constraint value at this point, and M_m is the number of points at which NO₃-N concentrations are monitored.

4.3 Risk

The degree to which constraints are respected depends on the decision-making stance with respect to risk. With a risk aversion value of 1.0 (see above), NO₃-N concentration constraints at all locations must be respected for all 475 realizations of posterior parameter fields employed by the northern Pineland Sands model. The number of realizations that can exceed any constraint increases in inverse proportion to the risk stance. Hence, for example, for a risk stance of 0.9, all constraints must be respected by 428 realizations out of the total of 475 realizations.

The lower the value of risk aversion, the greater the tolerance of NO₃-N concentrations in groundwater, and hence the greater the value of maximized NO₃-N leakage. As the latter is a surrogate for economic benefit, it follows that regional farm revenue grows with decreasing risk aversion. Risk in this study is an adjustable objective and maximized towards 1.0 (100% risk aversion), therefore serving as a competing objective against maximizing NO₃-N leakage.

4.4 Nitrate leaching rates and farm revenue

The optimization procedure which we have described thus far attempts to maximize regional farm revenue while respecting constraints on NO₃-N concentrations at certain points. Respect for these constraints leads to a reduction in the values of decision variables at certain locations. Recall that decision variables are NO₃-N leaching rates at 271 pilot points. These leaching rates are interpolated to farmed model cells. From Table 2, the leaching rates at these model cells for a given decision variable set are assigned a crop rotation based on the ranges each cell’s value falls within. The lower bound of the range represents the measured leaching rate for each crop rotation obtained from analyses undertaken by Errebhi et al. (1998), Wayment (2021), and Wilson et al. (2008a) as shown in Table 3. Subsequentially, a dollar yield per cell is obtained using the yield amount and economic value for each crop rotation listed in Table 3, with the economic value multiplied by the farmed area of each cell. This is then summed over all model cells to provide an estimate of regional farm revenue. Yields are converted to dollars per hectare using the 2023 average annual price per unit for each crop (dollars per bushel for corn and soybean, dollars per pound for dry bean, and dollars per hundredweight for potato) in Minnesota obtained from the Center for Farm FINBIN (2025), United States Department of AgricultureRisk Management Agency, 2024, United States Department of AgricultureNational Agriculture Statistics Service, 2024.

Table 2

Table 2. Optimal land use analysis upper and lower bounds for each crop rotation designation.

Table 3

Table 3. Crop rotations with associated NO3-N leaching rates and yield values (Errebhi et al., 1998; Wayment, 2021; Wilson et al., 2008a; Wilson et al., 2008b, FINBIN, 2025; United States Department of Agriculture Risk Management Agency, 2024; United States Department of Agriculture National Agriculture Statistics Service, 2024).

Obviously, this methodology for conversion of NO₃-N leaching rate to regional farm revenue is crude. Nevertheless, it satisfies the purpose of the present paper which is to demonstrate how modelling, combined with uncertainty analysis and constrained optimization, can inform a decision-making process while formulating a basis for negotiation between those who desire competing land-management objectives. Details of the modelling process, and the means by which objectives are computed, and constraints are enforced, can readily be varied to suit different decision contexts.

4.5 A social objective

As is demonstrated in the next section, respect for constraints imposed at NO₃-N monitoring points requires a lowering of NO₃-N leaching rates, and with it farm revenue. However, there may be more than one way in which constraints can be satisfied. In this demonstration study, we adopt as an additional management goal that the costs of respecting groundwater quality constraints be shared among many farmers to the extent that this is possible, rather than being borne by a few. Presumably, “the few” would be those whose farms are situated closest to monitoring points.

We therefore introduce another objective function for which a tradeoff can be sought with maximization of regional farm revenue on the one hand, and respect for monitoring constraints on the other hand. We name this the “farmer hardship penalty” and define it using Equation 5.

${\mathrm{\Phi }}_H={\displaystyle \sum _{i=1}^{M_d}}\mathit{max}\left(N_i-N_t,0\right)(5)$

In Equation 5 N_t is set to 61 kg NO₃-N ha^-1 yr^-1. This value corresponds to the NO₃-N leaching rate of corn-soybean crop rotation in Table 3. Below this leaching rate threshold, it is assumed that NO₃-N reducing practices need to be particularly stringent. Consequently, this is construed as the point at which farmers face economic hardship if NO₃-N targets are to be respected. The maximum value that Φ_H can adopt is 0.0. Under these circumstances it is conjectured that no farmer endures undue economic hardship.

5 Optimization outcomes

5.1 Results

Figure 7 shows Pareto fronts for NO₃-N concentration constraints of 5.4 mg L^-1 and 8 mg L^-1 at the monitoring wells of Figure 2. These were achieved after 74 iterations of the constrained optimization process at a cost of 62,625 model runs. Each dot in these figures corresponds to a set of decision variables; that is, each dot corresponds to a set of NO₃-N leaching rates ascribed to the 271 pilot points of Figure 2.

Figure 7

Figure 7. Pareto fronts calculated by PESTPP-MOU. Two Pareto fronts depicting the tradeoffs between the three objectives are shown, one for each scenario. Point colors denote each decision variable set’s farmer hardship penalty value, with darker colors reflecting lower values for this objective or increased farmer hardship. The decision variable sets marked by diamonds indicate those selected for further analysis.

It is immediately apparent from these plots that as the level of risk aversion rises, farm revenue falls. As farm revenue falls, the farmer hardship rises. Furthermore, the lower the NO₃-N non-exceedance threshold at monitoring points, the greater the penalty on farm revenue and the greater the economic hardship that must be endured to respect these thresholds.

Figures 7A,B suggest that at a risk level of 1.0, respect for monitoring thresholds cannot be guaranteed under any circumstances. This conclusion is obviously erroneous, as pervasive NO₃-N leaching rates of 0.0 would eliminate NO₃-N in groundwater and hence guarantee respect for all monitoring thresholds. These plots do not extend sufficiently far to the left to demonstrate this as an outcome of the finite number of decision variable points that are used to populate, and hence, define the Pareto front. Use of a greater number of points would overcome this problem, but at an intolerably high numerical cost.

Another feature of these plots that is somewhat anomalous is that the hardship index is negative (rather than zero) at a risk averse level of zero. There is little doubt that if the constrained optimization process were run for a few more iterations, then this would not be the case. However, this contradiction is tolerated because of the high numerical burden associated with defining the Pareto front.

Both plots in Figure 7 exhibit a relatively flat area to their left and a steeper zone to their right. This implies that there is a point of diminishing return whereby a risk-free guarantee of monitoring threshold adherence comes at a heavy economic cost and a high degree of farmer hardship. Conversely, a small to moderate degree of risk tolerance decreases farmer hardship dramatically.

The steepness of the right segments of all of these curves implies that dramatic mitigation of risks posed to groundwater quality can be implemented with little cost. However, as is demonstrated below, in spite of our efforts to distribute hardship through objective function definition, this cost is unduly borne by a comparatively small number of farms that are close to monitoring points. Once this “low hanging fruit” has been picked, reductions in NO₃-N leaching must be affected over larger areas by larger numbers of farmers to reduce monitored NO₃-N concentrations. With greater distance between farms and monitoring points, the effects of NO₃-N reduction at a particular farm on NO₃-N concentrations at a particular monitoring point become less certain. This requires more and more severity of NO₃-N leaching rate reduction and broadening of NO₃-N leaching rate reduction to gain surety that all monitoring thresholds are simultaneously respected.

To identify farming areas that are most impacted by the necessity to respect NO₃-N thresholds at monitoring points, we choose three points along the Pareto fronts. These points are depicted in Figure 7. Recall that each point in these graphs denotes a set of decision variables, that is, a set of NO₃-N leaching rates ascribed to pilot points. These leaching rates are spatially interpolated to model cells. Leaching rates pertaining to three levels of risk for the NO₃-N concentration constraints (shown in Figure 2) are presented in Figure 8. Cooler colors denote greater perturbations from a default NO₃-N leaching rate of 150 kg N ha^-1 and hence greater reductions in farm revenue. Included in each of these figures is a pie chart. This denotes the subdivision of present-day farming areas into different cropping types based on optimized NO₃-N leaching rate; the contents of Table 2 are used to convert NO₃-N leaching rate to crop type.

Figure 8

Figure 8. NO₃ -N leaching rate optimization results for the selected solutions accompanied by their optimal land use by percent area pie charts. (A–C) corresponds to the 5.4 mg L^-1 constraint scenario while (D–F) correspond to the 8 mg L^-1 constraint scenario The percentage area of the primary crops grown in the region for 2024 were 13% potato, 30.8% corn, 27.1% dry bean, 12.7% soybean. The remaining 16.3% of the total cropland was other minor crop types such as wheat, oats, alfalfa, peas, and sugarbeet (United States Department of Agriculture, National Agriculture Statistics Service, 2025).

Figure 8 makes it clear that the lower the threshold of NO₃-N concentration acceptability at monitoring points, the greater the need for alteration of farming practices from the highest revenue generating rotation to comply with these thresholds. Areas for which greatest adjustments are required tend to be those that are closest to NO₃-N monitoring points. Of particular note are two central regions of cultivated cropland, near the Straight River, City of Park Rapids, and the Shell River. Additionally, a small patch of lower optimal leaching rates appears in the northwest section of the model domain within the Pine Point Township on White Earth Reservation.

Table 4 attempts to provide some economic context to the maps of Figure 8. This analysis must be considered illustrative at best. However, it suggests that the economic cost of safeguarding the northern Pineland Sands Aquifer System is far from trivial. It is conjectured that the same applies to many other aquifer systems whose NO₃-N levels are high as a result of agricultural activities.

Table 4

Table 4. Optimal solution land use and economic analysis results for each constraint scenario.

5.2 Stakeholder preferences

Results from multi-objective optimization (MOO) are preference-agnostic as the Pareto front makes no inference on which objectives matter more to decision-makers. On the contrary, stakeholders typically have preferences for outcomes which favor objectives that matter the most to them. Multi-criteria analysis (MCA) can be used in the post-processing of Pareto fronts to rank each optimal solution in accordance with the preferences of each stakeholder group. MCA has been previously applied in conjunction with MOO in the water resources field (Almasri and Kaluarachchi, 2005; Alizadeh et al., 2017; Gaur et al., 2021).

In this framework, the MCA ranking method, TOPSIS, is applied to the non-dominated solutions resulting from the 8 mg L^-1 constraint scenario. Each objective is vector-normalized in accordance with the guidance on TOPSIS outlined in Katanski (2017). TOPSIS operates under the idea that the best solution out of a set of alternatives will be closest to the “ideal best solution” and furthest away from the “ideal worst solution”.

In this demonstration, three stakeholder groups are assumed, representing local growers, environmental stakeholders, and local government agencies. Hypothetical weights are then assigned to each objective based on assumed stakeholder preferences being:

• Local growers care three times as much about the economic objective (average decision variable leaching rate), and twice as much about the social objective (farmer hardship), then the environmental objective (risk aversion).

• Environmental stakeholders care three times more about the environmental objective than the social or economic objectives

• Local government cares half as much about the economic objective than the environmental and social objectives

These preferences are translated into weights (Table 5), where for each stakeholder group, the sum of weights must equal 1.

Table 5

Table 5. Objective weights derived from hypothetical stakeholder preferences.

For each stakeholder group, TOPSIS then ranks each solution based on the weighted normalized objective values using Equations 6a–6c.

$S_i^{*}=\sqrt{{\displaystyle \sum _{j=1}^J}{\left(v_{ij}-v_j^{*}\right)}^2}(\mathrm{6a})$

$S_i^{-}=\sqrt{{\displaystyle \sum _{j=1}^J}{\left(v_{ij}-v_j^{-}\right)}^2}(\mathrm{6b})$

$V_i=\frac{S_i^{-}}{S_i^{*}+S_i^{-}}(\mathrm{6c})$

Where $S_i^{*}$ is the Euclidean distance for alternative i to the ideal best solution, $S_i^{-}$ is the Euclidean distance for alternative i to the ideal worst solution, $v_{ij}$ is the weighted, normalized values of alternative i with respect to object j, $v_j^{-}$ is the ideal worst value for objective j, $v_j^{*}$ is the ideal best value for objective j, and $V_i$ is the relative distance to the ideal point for each alternative solution.

TOPSIS results in rankings for each alternative based on a stakeholder’s preference. When multiple stakeholders are involved, these rankings can vary depending on what matters most to each group. Borda score voting (Borda, 1781; Black, 1958) is a social choice rule methodology which determines a winning solution based on aggregated voter rankings for a given set of alternative solutions. This methodology assigns a score based on the premise in Equation 7:

$S_{ij}={\displaystyle \sum _{v=1}^V}n-r_{ij}+1(7)$

Where S_ij is the Borda score for an alternative solution, n is the number of solutions, and r_ij is a solutions rank, and V is the number of voters which in this context is stakeholder groups. Applying Borda score voting to the ranked alternative solutions from resulting from the previous analysis, the winning solution given the hypothetical preferences of each stakeholder group has as average decision variable value of 138.67 kg NO₃-N ha^-1 yr^-1, a farmer hardship value of −674, and a risk aversion of 0.789. The resulting optimal maximum leaching rate field is shown in Figure 9.

Figure 9

Figure 9. The Borda count winning solution based on the hypothetical stakeholder preferences assigned previously to each objective. This solution has a risk aversion of 0.789, a farmer hardship index of −674, and a average decision variable value of 138.67 kg NO₃-N ha^-1 yr^-1.

As one can observe, this information pipeline demonstrated within can be readily adapted to include more stakeholder groups, and in doing so, may help further increase stakeholder engagement in the decision-analysis process by allowing a practitioner to factor in the preferences of those most impacted by a management outcome in their attempt at finding the “least objectionable decision” for a given application of this framework.

6 Discussion

6.1 Limitations

We begin our discussion of analyses that are presented in the previous section by listing a few obvious limitations. The intent of this paper is primarily to demonstrate a methodology for processing environmental data in a manner that can add value to discussions about land management in general, and management of the northern Pineland Sands Aquifer System in particular. Outcomes of our analysis should be taken as illustrative rather than definitive. (We also note that no methodology for predicting the future behavior of a groundwater system, and for quantifying the uncertainty of that behavior, is immune from a plethora of limitations).

The groundwater model on which our analyses are based is steady state. This has two repercussions. The first is that its predictions take no account of variation of NO₃-N concentration with time. The second is that accommodation of time varying groundwater levels in history-matching of a steady state model requires recognition of “measurement noise” that accounts for errors in the temporal averaging process. This can inflate the uncertainties of groundwater model parameters, together with predictions that are sensitive to these parameters. History-matching and prediction using a transient model would, however, have added considerably to model run times. Furthermore, it is not a foregone conclusion that parameter and predictive uncertainties would have been greatly reduced through transient history-matching. This is because transient history-matching must inform many more parameters than steady state history-matching.

The simulation process faces other shortcomings as well. Uncertainties in groundwater denitrification rates in groundwater have been considered in only a cursory manner. As has been discussed, a lower-than-assumed groundwater denitrification rate would result in greater alterations to management practices than is indicated herein. Greater denitrification rates would do the opposite. As is explained above, we adopted a conservative option. It is worth remarking, however, that high NO₃-N concentrations found in historic NO₃-N data for the region are commensurate with model predictions that are presented herein.

As discussed above, optimization under uncertainty requires that many model runs be undertaken, especially where these runs are devoted to definition of a Pareto front. Outcomes of such an analysis may depend to some extent on the number of model runs that are undertaken, and hence on the number of iterations over which the optimization process is allowed to proceed. There are indications in the results that are presented herein that Pareto fronts may have been shifted slightly upwards and outwards after more optimization iterations. However, it is most unlikely that the shapes of these fronts would have changed, nor the locations and magnitudes of alterations to land use that are identified as being required to satisfy environmental monitoring constraints.

Most importantly, the link between NO₃-N leakage rate and farm revenue that we adopt in our analysis is approximate at best. Hence the economic aspects of our analysis are inexact. Adoption of this crude linkage allowed formulation of a relatively simple cost function. It is important to note, however, that in contrast to simulation of groundwater and contaminant movement, calculation of cost does not, in general, incur a high computational burden. Hence the use of alternative costs functions in constrained optimization-under-uncertainty is not a difficult matter. Alternative cost functions may include the introduction of measures to reduce NO₃-N leaching that do not require alteration to crop type. The benefit of using a simple cost function in analyses presented above is that it readily identifies locations where the need for leaching mitigation practices is greatest, and where the need to observe environmental constraints may render investments in these practices worthwhile.

6.2 Remarks on methodology

As noted in the introduction, while groundwater modelling has the potential to make significant contributions to scientifically based environmental management, it has not always fulfilled expectations. Ideally, introduction of the scientific method to discussions between parties who hold different values (and hence view optimality of environmental management through different lenses) should create a space for informed and objective discussion that would not otherwise take place. The scientific method demands that those who visit this space attempt to explore the repercussions of different land management options based on facts and data alone. Just as importantly, it demands that consequences of data insufficiency be fully accommodated in all analyses, and that prognostications of future system behavior be therefore accompanied by uncertainty intervals. Scientifically-based data processing should attempt to reduce these intervals to the extent that available data allows, but no further.

While decision-stakeholders may hold differing opinions on matters such as the nature of the prior parameter probability distribution, and some aspects of model structure, in principle these matters are resolvable. Where they cannot be resolved, they can be accommodated through adjustments to prior probability distributions that meet these concerns. Furthermore, though not discussed herein, these probability distributions can be categorical in nature if this is required, therefore allowing for discontinuous possibilities for some aspects of system behavior, and/or disputed or unknown subsurface properties.

Stakeholder values are introduced to a decision-making process during a formal or informal optimization process. In the study that is documented herein, optimization is formal. In this case, as in all other cases, the manner in which the optimization problem is formulated, and the constraints that are imposed as an objective function is maximized or minimized, are designed to reflect stakeholder values and concerns. Using the optimization-under-uncertainty methodology that is demonstrated herein, the repercussions of proposed management options on valued aspects of system behavior can be readily examined. So too can the imposition of monitoring thresholds, and the respect that is given to these thresholds according to different proposed system management protocols. The utility of this methodology is enhanced by its ability to define tradeoff curves which delineate the stakeholder-perceived costs and benefits of pursuing different management options. Identification of points of diminishing return along these tradeoff curves may have significant impacts on stakeholder negotiations.

Groundwater modelling on its own may assist stakeholders in understanding system behavior to some extent. However, value must be added to the modelling process for it to support the making of decisions, particularly in contexts where no management option is likely to please all parties. This study illustrates two important numerical activities, with simulation at their heart, that can provide considerable support to the decision-making process. These activities do not absolve a decision-maker from the onerous task of having to choose between different options that may be valued differently by different stakeholder groups. However, when undertaken together, they provide an intellectual framework which enables separation of science from values. This protects the scientific process from criticism, at the same time as it preserves its ability to create a unique and important space in which discussions can be held between parties for whom such discussions may otherwise be impossible. Should the decision-making process dissolve into legal wrangling, the importance of separating science from values becomes even more important.

The first of these two activities is simulator-based environmental data processing. In the study that is documented herein, this is implemented through Bayesian history-matching of observed groundwater behavior. Through this process, information is assimilated by a model’s parameters. This information can then be transferred to model predictions as management options are tested. If properly implemented, the outcome of simulator-based data assimilation is an ability to make bias-free predictions of system behavior under a variety of different management conditions. Importantly, these predictions are accompanied by estimates of their uncertainties. In some management contexts, these predictions may involve direct comparisons of the outcomes of different management options. Experience and anecdotal evidence repeatedly reveal that uncertainties associated with predictive differences are smaller than those associated with direct predictions of system states and fluxes. In either case, quantification of these uncertainties is essential to definition and management of the risks associated with any contemplated course of management action.

The second of these activities is decision optimization. Conceptually, an optimization process can be formulated in a way that meets the concerns of any stakeholder group. This is achieved through definition of appropriate objective functions, through formulation of constraints that are imposed on these functions, and through definition of the severity with which these constraints are imposed. The methodology that is presented herein allows tradeoffs between the interests of different stakeholder groups to be explored. Importantly, one of these tradeoffs is acceptance or aversion to risk.

As previously discussed, decision-making is rarely value-free. Instead, most stakeholders hold specific values which may give them a predisposition towards one outcome or another based on which objective said outcome favors the most. While the two prior activities allow for the separation of science from values, the third activity shown assists in bringing the values back into the decision-making process. The pairing of MCA with the methodology shown above facilitates an information pipeline which can aid in the difficult task of finding the “compromise” out of the many alternatives found during the optimization process. While most likely no management decision will satisfy all parties, some alternatives may be less objectionable to all parties than others. The inclusion of MCA into this framework, to the authors’ belief, can help enable further stakeholder engagement by factoring in the preferences of those most affected by a management outcome into the decision-making process.

For the northern Pinelands Sands Aquifer System, the study that is document herein demonstrates that to achieve a high degree of risk aversion for the NO₃-N concentration constraint locations and thresholds used in this analysis, “business as usual” is unlikely to be an acceptable management option. Unfortunately, it also suggests that stakeholders may hold differing opinions on alternative management options that are proposed to ameliorate risks to human and biotic health that have emerged from current farming practices. These options may differ in their identification of.

• The economic burden that must be borne;

• The distribution of that burden;

• The risks that remain to waterways and human health after alterations to farming practices have been made.

The study suggests that reductions of farmland from the higher-impact, higher-grossing crop rotations may be required to respect the constraints used in this analysis. It also shows where these alterations are likely to be most effective in improving groundwater quality. However, it must be borne in mind that the spread of colors in maps such as Figure 8 is strongly affected by the disposition of NO₃-N concentration constraints. With the introduction of more constraints to the optimization process, requirements for reductions on NO₃-N leaching may be more widespread than is indicated in Figure 8. Decisions on where extra monitoring points should be located, thresholds that are ascribed to these monitoring points, and the severity with which these thresholds should be respected are all matters that must be informed by stakeholder values. Management options and tradeoffs that arise from these values can be readily explored using methodologies such as those that are demonstrated herein.

7 Conclusion

7.1 Summary of findings

Provisional, simulator-based exploration of the trade-off between regional farm production and achievement of water quality objectives in the Pineland Sand Aquifer System has demonstrated that alterations to farming practices from higher-impact crop rotations may be necessary if NO₃-N water quality objectives are to be met. However, there is some degree of uncertainty associated with the nature and extent of these alterations. This is an outcome of (among other things) uncertainties in aquifer hydraulic properties and in regional denitrification rates. Enforcement of changes to farming practices therefore requires that decision-makers accept a level of risk that these changes will not achieve desired objectives. This level of risk is incorporated in Pareto tradeoff curves that emerged from our study.

Nevertheless, our analyses have demonstrated that significant improvements in groundwater quality at key monitoring points can be achieved with a comparatively small decrease in agricultural yield in the general study area. However, in spite of an equity goal introduced to our analysis that the burden of change be spread among many farmers rather than resting on just a few, it is apparent that greatest gains in water quality can be achieved by reducing NO₃-N leaching rates close to monitoring points.

7.2 Implications for science-based decision-making

Analyses that are described herein support participatory decision-making in a number of ways. First, numerical modelling is deployed in a probabilistic way that encapsulates as much hydrogeological knowledge of a study area as is currently available. Limitations of current knowledge, and consequences for predictive uncertainty, are therefore recognized. Second, model-based history-matching to informative data can reduce uncertainties of decision-critical predictions through assimilation of information that is contained in measurements of historical system behavior. All of this is done in ways that are in accordance with science-based understanding of aquifer properties and groundwater behavior in the study area.

Nevertheless, while maintaining its scientific basis, methodologies that are described herein are welcoming of input by different groups with disparate interests. First, stakeholders are able to contribute to the knowledge base on which the numerical model is built; where parties disagree in decision-salient hydrogeological concepts, these disagreements can be expressed in model parameterization uncertainties that are integral to model deployment. Second, the issues that modelling addresses should reflect stakeholder concerns. In the present study, the focus is on water quality at specified monitoring points. However other issues could be easily addressed with the model described herein, or with other models that are dedicated to those issues. Third, analyses described herein make no assumptions on how the interests of one stakeholder group should be balanced against those of another. Examination of the repercussions of a broad range of management strategies is built into these analyses. Fourth, the analyses recognize that outcomes of any alteration to regional management are uncertain, and that managers must acknowledge the risk that any decision that they take may not have the outcomes that they desire. It is interested parties, and not disinterested scientists, that decide on the level of risk that will determine decisions that are ultimately made. Meanwhile, the repercussions of these risks are illuminated by scientific analysis. Fifth, the preferences of each stakeholder can be considered numerically and utilized to rank and vote on each alternative management outcome to attempt to find a “middle-ground” between all affected parties. The “winning” solution is dependent on the weights that decision-makers assign to each objective, therefore, leaving the ranking of alternatives to the stakeholders rather than the scientists who conduct this analysis.

Decisions are taken by people. Rarely do people yield their authority as stakeholders or as decision-makers to an optimization algorithm, or to a tradeoff curve, regardless of the objectivity on which tradeoff calculations are based. At the same time, stakeholders are eager to claim that their decisions are based on science. The interplay between science and decision-making is therefore complex.

While it is not the purpose of this paper to explore this interplay, it is hoped that methodologies described herein have the ability to exert a positive impact on this interplay.

7.3 Future applications and limitations

Analyses that are discussed herein have general applicability. They have been conducted with public domain software using well-documented algorithms. They are therefore readily applied to other sites and to other issues.

Ideally, future applications of this methodology in the Pineland Sands Aquifer should be accompanied by community consultation. Part of the consultation process would be a demonstration of what analyses such as those that are presented herein are able to contribute to participatory decision-making. Likewise, the weighting of each objective relative to the others can only be accomplished through engagement with stakeholders most effected by the outcomes of this methodology described herein. The next steps are up to decision participants. However, given some of the findings that have emerged from the demonstration study that is presented herein, it is possible to make some suggestions.

It has been demonstrated that, for many of the monitoring points at which NO₃-N concentration limits may be imposed, those whose farms are close to the points exert the largest impact. Management of the regional aquifer may therefore be better served by a series of local models rather than by a single, large model. Development of software that enables rapid construction, deployment and history-matching of these models therefore becomes important. At the same time, these models may be served by more realistic representation of alluvial hydraulic property connectedness than was used in present modeling. This can be achieved using recently developed methods that incorporate non-stationary geostatistics; see Kitlasten et al. (2005).

The extent to which methodologies that are described herein can impact discussions between different stakeholder groups, and the ability of these groups to participate in decisions which gain the commitment of all parties, is yet to be seen. However, it is hoped that their sound basis in science, together with their porous interface to stakeholder input, will grant them a positive impact. Documentation of this impact awaits future publications.

Data availability statement

The datasets presented in this article are not readily available because the locations of some constraints used in this study are restricted as per Minnesota Department of Health policy. To obtain the data used in this study, fill out the form at the following link and contact the corresponding author upon receiving approval. Requests to access the datasets should be directed to https://survey.vovici.com/se/56206EE33CCF9BC4 and JN, bmllYmVyQHVtbi5lZHU=.

Author contributions

PM: Investigation, Writing – review and editing, Conceptualization, Software, Validation, Writing – original draft, Data curation, Visualization, Project administration, Methodology, Formal Analysis. JD: Software, Writing – review and editing, Methodology, Validation, Writing – original draft, Conceptualization, Supervision. LL: Visualization, Writing – review and editing, Data curation. JN: Supervision, Funding acquisition, Resources, Writing – review and editing. BW: Supervision, Writing – review and editing, Conceptualization, Methodology. JM: Resources, Writing – review and editing, Funding acquisition, Supervision.

Funding

The authors declare that financial support was received for the research and/or publication of this article. This study was financially supported by the Legislative-Citizen Commission on Minnesota Resources (M.L. 2023, chp 60, art. 2, sec. 2 subd. 04m).

Acknowledgements

We would like to thank our partners at the White Earth Band of Minnesota Chippewa Tribe. We would also like to acknowledge George Kraft from the University of Wisconsin-Stevens Point, Anthony Runkle from the Minnesota Geological Survey, and Eduardo Garay Lagos from the University of Minnesota for their guidance on this study.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The authors declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenvs.2025.1634272/full#supplementary-material

References

Alizadeh, M. R., Nikoo, M. R., and Rakhshandehroo, G. R. (2017). Developing a multi-objective conflict-resolution model for optimal groundwater management based on fallback bargaining models and social choice rules: a case study. Water Resour. Manag. 31 (5), 1457–1472. doi:10.1007/s11269-017-1588-7

CrossRef Full Text | Google Scholar

Almasri, M. N., and Kaluarachchi, J. J. (2005). Multi-criteria decision analysis for the optimal management of nitrate contamination of aquifers. J. Environ. Manag. 74 (4), 365–381. doi:10.1016/j.jenvman.2004.10.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Almasri, M. N., and Kaluarachchi, J. J. (2007). Modeling nitrate contamination of groundwater in agricultural watersheds. J. Hydrology 343 (3-4), 211–229. doi:10.1016/j.jhydrol.2007.06.016

CrossRef Full Text | Google Scholar

Ameur, M., Aouiti, S., Hamzaoui-Azaza, F., Cheikha, L. B., and Gueddari, M. (2021). Vulnerability assessment, transport modeling and simulation of nitrate in groundwater using SI method and modflow-MT3DMS software: case of sminja aquifer, Tunisia. Environ. Earth Sci. 80, 220. doi:10.1007/s12665-021-09491-z

CrossRef Full Text | Google Scholar

Arbuckle, T. E., Sherman, G. J., Corey, P. N., Walters, D., and Lo, B. (1988). Water nitrates and CNS birth defects: a population-based case-control study. Archives Environ. Health An Int. J. 43 (2), 162–167. doi:10.1080/00039896.1988.9935846

PubMed Abstract | CrossRef Full Text | Google Scholar

Aschebrook-Kilfoy, B., Heltshe, S. L., Nuckols, J. R., Sabra, M. M., Shuldiner, A. R., Mitchell, B. D., et al. (2012). Modeled nitrate levels in well water supplies and prevalence of abnormal thyroid conditions among the old order amish in Pennsylvania. Environ. Health 11, 6. doi:10.1186/1476-069X-11-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Bauer, E., Chandler, V. W., Jirsa, M., Marshall, K., Gowan, A., Hamilton, J., et al. (2016). Data from: C-42, geologic atlas of becker county, Minnesota. Paul, MN: University Digital Conservancy. Available online at: https://hdl.handle.net/11299/184650.

Google Scholar

Bayer, P., de Paly, M., and Bürger, C. M. (2010). Optimization of high-reliability-based hydrological design problems by robust automatic sampling of critical model realizations. Water Resour. Res. 46 (5), W05504. doi:10.1029/2009WR008081

CrossRef Full Text | Google Scholar

Bedekar, V., Morway, E. D., Langevin, C. D., and Tonkin, M. J. (2016). MT3D-USGS version 1: a US geological survey release of MT3DMS updated with new and expanded transport capabilities for use with MODFLOW. U.S. Geol. Surv. Tech. Methods, 6–A53. doi:10.3133/tm6A53

CrossRef Full Text | Google Scholar

Black, D. (1958). The theory of committees and elections. Cambridge: Cambridge University Press.

Google Scholar

Borda, J. D. (1781). M'emoire sur les' elections au scrutin. Hist. l'Acad'emie R. Sci. 102, 657–665.

Google Scholar

Bradt, R. J. (2023). “Groundwater atlas of becker county, Minnesota: part B hydrogeology,” in County atlas series C-42. St. Paul, MN: Minnesota Department of Natural Resources, 3.

Google Scholar

Brender, J. D., Olive, J. M., Felkner, M., Suarez, L., Marckwardt, W., and Hendricks, K. A. (2004). Dietary nitrites and nitrates, nitrosatable drugs, and neural tube defects. Epidemiology 15 (3), 330–336. doi:10.1097/01.ede.0000121381.79831.7b

PubMed Abstract | CrossRef Full Text | Google Scholar

Brender, J. D., Weyer, P. J., Romitti, P. A., Mohanty, B. P., Shinde, M. U., Vuong, A. M., et al. (2013). Prenatal nitrate intake from drinking water and selected birth defects in offspring of participants in the national birth defects prevention study. Environ. Health Perspect. 121 (9), 1083–1089. doi:10.1289/ehp.1206249

PubMed Abstract | CrossRef Full Text | Google Scholar

Budde, N. R. (2024). “Groundwater atlas of hubbard county, Minnesota: part B hydrogeology,” in County atlas series C-41. St. Paul, MN: Minnesota Department of Natural Resources, 3–pls.

Google Scholar

Budde, N. R., and Baratta-Person, V. M. (2024). “Groundwater atlas of wadena county, Minnesota: part B hydrogeology,” in County atlas series C-40. 2. St. Paul, MN: Minnesota Department of Natural Resources.

Google Scholar

Cartwright, N., Clark, L., and Bird, P. (1991). The impact of agriculture on water quality. Outlook Agric. 20 (3), 145–152. doi:10.1177/003072709102000304

CrossRef Full Text | Google Scholar

Chen, Y., and Oliver, D. S. (2013). Levenberg marquardt forms of the iterative ensemble smoother for efficient history matching and uncertainty quantification. Comput. Geosci. 17 (4), 689–703. doi:10.1007/s10596-013-9351-5

CrossRef Full Text | Google Scholar

Correa-González, A., Hernández-Bedolla, J., Martínez-Cinco, M. A., Sánchez-Quispe, S. T., and Hernández-Hernández, M. A. (2023). Assessment of nitrate in groundwater from diffuse sources considering spatiotemporal patterns of hydrological systems using a coupled SWAT/MODFLOW/MT3DMS model. Hydrology 10 (11), 209. doi:10.3390/hydrology10110209

CrossRef Full Text | Google Scholar

Coulon, C., Lemieux, J.-M., Pyret, A., Bayer, P., Young, N. L., and Molson, J. (2022). Pumping optimization under uncertainty in an island freshwater lens using a sharp-interface seawater intrusion model. Water Resour. Res. 58 (8), e2021WR031793. doi:10.1029/2021WR031793

CrossRef Full Text | Google Scholar

Coulon, C., White, J. T., Pryet, A., Gatel, L., and Lemieux, J.-M. (2024). An ensemble-based approach for pumping optimization in an island aquifer considering parameter, observation and climate uncertainty. Hydrol. Earth Syst. Sci. 28, 303–319. doi:10.5194/hess-28-303-2024

CrossRef Full Text | Google Scholar

Croen, L. A., Todoroff, K., and Shaw, G. M. (2001). Maternal exposure to nitrate from drinking water and diet and risk for neural tube defects. Am. J. Epidemiol. 153 (4), 325–331. doi:10.1093/aje/153.4.325

PubMed Abstract | CrossRef Full Text | Google Scholar

Daniel, T. J., Richendrfer, J., Falta, R., Murdoch, L., Lin, H., and Darnault, C. J. (2023). Numerical modeling of groundwater flow and nitrate transport in karst and wastewater irrigated agricultural and forest landscapes, PA, USA. Agrosystems, Geosciences and Environ. 6, e20427. doi:10.1002/agg2.20427

CrossRef Full Text | Google Scholar

De Roos, A. J., Ward, M. H., Lynch, C. F., and Cantor, K. P. (2003). Nitrate in public water supplies and the risk of colon and rectum cancers. Epidemiology 14 (6), 640–649. doi:10.1097/01.ede.0000091605.01334.d3

PubMed Abstract | CrossRef Full Text | Google Scholar

Deb, K., Pratap, A., Agarwal, S., and Meyarivan, T. A. M. T. (2002). A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Trans. Evol. Comput. 6 (2), 182–197. doi:10.1109/4235.996017

CrossRef Full Text | Google Scholar

Dewitz, J. (2021). Data from: National land cover database (NLCD) 2019 products (ver. 2.0, June 2021). U.S. Geol. Surv. Data Release. doi:10.5066/P9KZCM54

CrossRef Full Text | Google Scholar

Doherty, J. (2015). Calibration and uncertainty analysis for complex environmental models. Brisbane, Australia: Watermark Numerical Computing.

Google Scholar

Doherty, J., and Simmons, C. (2013). Groundwater modelling in decision support: reflections on a unified conceptual framework. Hydrogeology J. 21, 1531–1537. doi:10.1007/s10040-013-1027-7

CrossRef Full Text | Google Scholar

Dorsch, M. M., Scragg, R. K., McMichael, A. J., Baghurst, P. A., and Dyer, K. F. (1984). Congenital malformations and maternal drinking water supply in rural South Australia: a case-control study. Am. J. Epidemiol. 119 (4), 473–486. doi:10.1093/oxfordjournals.aje.a113764

PubMed Abstract | CrossRef Full Text | Google Scholar

Eberhart, R., and Kennedy, J. (1995). “A new optimizer using particle swarm theory,” in Proceedings of the sixth international symposium on micro machine and human science (Nagoya, Japan), 39–43. doi:10.1109/MHS.1995.494215

CrossRef Full Text | Google Scholar

Errebhi, M., Rosen, C. J., Gupta, S. C., and Birong, D. E. (1998). Potato yield response and nitrate leaching as influenced by nitrogen management. Agron. J. 90 (1), 10–15. doi:10.2134/agronj1998.00021962009000010003x

CrossRef Full Text | Google Scholar

Espejo-Herrera, N., Gràcia-Lavedan, E., Boldo, E., Aragonés, N., Pérez-Gómez, B., Pollán, M., et al. (2016). Colorectal cancer risk and nitrate exposure through drinking water and diet. Int. J. Cancer 139 (2), 334–346. doi:10.1002/ijc.30083

PubMed Abstract | CrossRef Full Text | Google Scholar

Fienen, M. N., Haserodt, M. J., Leaf, A. T., and Westenbroek, S. M. (2022). Simulation of regional groundwater flow and groundwater/lake interactions in the central sands. Wisconsin: U.S. Geological Survey Scientific Investigations Report 2022-5046. doi:10.3133/sir20225046

CrossRef Full Text | Google Scholar

Fienen, M. N., Corson-Dosch, N. T., Jahn, K. L., and White, J. T. (2024). Comparing single and multiple objective constrained optimization algorithms for tuning a groundwater remediation system. Environ. Model. and Softw. 173, 105952. doi:10.1016/j.envsoft.2024.105952

CrossRef Full Text | Google Scholar

FINBIN (2025). “Center for farm financial management,”. Paul, MN: University of Minnesota. Available online at: http://finbin.umn.edu (Accessed April 24, 2025).

Google Scholar

Frind, E. O., Duynisveld, W. H. M., Strebel, O., and Böttcher, J. (1990). Modeling of multicomponent transport with microbial transformation in groundwater: the fuhrberg case. Water Resour. Res. 26 (8), 1707–1719. doi:10.1029/WR026i008p01707

CrossRef Full Text | Google Scholar

Gatseva, P. D., and Argirova, M. D. (2008). High-nitrate levels in drinking water May be a risk factor for thyroid dysfunction in children and pregnant women living in rural Bulgarian areas. Int. J. Hyg. Environ. Health 211 (5-6), 555–559. doi:10.1016/j.ijheh.2007.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaur, S., Srinivasa Raju, K., Nagesh Kumar, D., and Bajpai, M. (2021). Multicriterion decision making in groundwater planning. J. Hydroinformatics 23 (3), 627–638. doi:10.2166/hydro.2021.122

CrossRef Full Text | Google Scholar

Gorelick, S. M. (1997). “Incorporating uncertainty into aquifer management models,” in Subsurface flow and transport. Editors G. Dagan, and S. P. Neuman (Cambridge University Press), 101–112.

CrossRef Full Text | Google Scholar

Harbaugh, A. W. (2005). MODFLOW-2005: the U.S. geological survey modular groundwater model-the groundwater flow process. U.S. Geol. Surv. Tech. Methods, 6–A16. doi:10.3133/tm6A16

CrossRef Full Text | Google Scholar

Katanski, A. V. (2017). Stories that nourish: minnesota anishinaabe wild rice narratives. Am. Indian Cult. Res. J. 41 (3). doi:10.17953/aicrj.41.3.katanski

CrossRef Full Text | Google Scholar

Kay, J., and King, M. (2020). Radical uncertainty: decision-making beyond the numbers. WW Norton and Company.

Google Scholar

Kincheloe, J. W., Wedemeyer, G. A., and Koch, D. L. (1979). Tolerance of developing salmonid eggs and fry to nitrate exposure. Bull. Environ. Contam. Toxicol. 23 (4-5), 575–578. doi:10.1007/BF01770006

PubMed Abstract | CrossRef Full Text | Google Scholar

Kitlasten, W., Moore, C., and Doherty, J. (2025). Exploring the use of new data assimilation technologies to map groundwater quality vulnerability in a large alluvial aquifer. Front. Earth. Sci. Sec. Hydrosphere. 13, 1609778. doi:10.3389/feart.2025.1609778

CrossRef Full Text | Google Scholar

Langevin, C. D., Hughes, J. D., Banta, E. R., Niswonger, R. G., Panday, S., and Provost, A. M. (2017). Documentation for the MODFLOW 6 groundwater flow model. U.S. Geol. Surv. Tech. Methods, 6–A55. doi:10.3133/tm6A55

CrossRef Full Text | Google Scholar

Llopis-Albert, C., Palacios-Marqués, D., and Merigó, J. M. (2014). A coupled stochastic inverse-management framework for dealing with nonpoint agriculture pollution under groundwater parameter uncertainty. J. Hydrology 511, 10–16. doi:10.1016/j.jhydrol.2014.01.021

CrossRef Full Text | Google Scholar

Lusardi, B. (2016). Data from: C-40, geologic atlas of wadena county. Minnesota: Paul, MN: University Digital Conservancy. Available online at: https://hdl.handle.net/11299/183206.

Google Scholar

Lusardi, B. (2018). Data from: C-41, geologic atlas of hubbard county, Minnesota. Paul, MN: University Digital Conservancy. Available online at: https://hdl.handle.net/11299/198898.

Google Scholar

Lyne, V., and Hollick, M. (1979). “Stochastic time-variable rainfall-runoff modelling,” in Paper presented at: Institute of engineers Australia national conference on hydrology and water resources. Perth, Australia, 89–93.

Google Scholar

McElroy, J. A., Trentham-Dietz, A., Gangnon, R. E., Hampton, J. M., Bersch, A. J., Kanarek, M. S., et al. (2008). Nitrogen-nitrate exposure from drinking water and colorectal cancer risk for rural women in Wisconsin, USA. J. Water Health 6 (3), 399–409. doi:10.2166/wh.2008.048

PubMed Abstract | CrossRef Full Text | Google Scholar

McGurk, M. D., Landry, F., Tang, A., and Hanks, C. C. (2006). Acute and chronic toxicity of nitrate to early life stages of Lake trout (Salvelinus namaycush) and Lake whitefish (Coregonus clupeaformis). An Int. J. 25 (8), 2187–2196. doi:10.1897/05-270r.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Minnesota Department of Agriculture (MDA) (2019). Groundwater protection rule, Minnesota administrative rules. Available online at: https://www.revisor.mn.gov/rules/1573/(Accessed April 23, 2025).

Google Scholar

Minnesota Department of Health (MDH) (2025). Data from: minnesota well index. Available online at: https://mnwellindex.web.health.state.mn.us/.

Google Scholar

Minnesota Department of Natural Resources (MNDNR) (2025). Data from: cooperative groundwater monitoring program. Available online at: https://www.dnr.state.mn.us/waters/cgm/index.html.

Google Scholar

Minnesota Department of Natural Resources (MNDNR) and U.S. Geological Survey (USGS) (2025). Data from: cooperative stream gaging. Available online at: https://www.dnr.state.mn.us/waters/csg/index.html.

Google Scholar

Minnesota Pollution Control Agency (MPCA) (2014). Minnesota nutrient reduction strategy. Available online at: https://www.pca.state.mn.us/sites/default/files/wq-s1-80.pdf (Accessed October 18, 2024).

Google Scholar

Monson, P. (2022). Aquatic life water quality standards draft technical support document for nitrate. Minn. Pollut. Control Agency. Available online at: https://www.pca.state.mn.us/sites/default/files/wq-s6-13.pdf (Accessed October 23, 2024).

Google Scholar

Niswonger, R. G., Panday, S., and Ibaraki, M. (2011). “MODFLOW-NWT, a newton formulation for MODFLOW-2005,” in U.S. geological survey, techniques and methods 6-A37. doi:10.3133/tm6A37

CrossRef Full Text | Google Scholar

Novotny, V. (1999). Diffuse pollution from agriculture-a worldwide outlook. Water Sci. Technol. 39 (3), 1–13. doi:10.2166/wst.1999.0124

CrossRef Full Text | Google Scholar

Oliver, D. S. (2022). Hybrid iterative ensemble smoother for history matching of hierarchical models. Math. Geosci. 54 (8), 1289–1313. doi:10.1007/s11004-022-10014-0

CrossRef Full Text | Google Scholar

Panday, S., Langevin, C. D., Niswonger, R. G., Ibaraki, M., and Hughes, J. D. (2013). MODFLOW–USG version 1: an unstructured grid version of MODFLOW for simulating groundwater flow and tightly coupled processes using a control volume finite-difference formulation. U.S. Geol. Surv. Tech. Methods 6-A45. doi:10.3133/tm6A45

CrossRef Full Text | Google Scholar

PEST Development Team (2025). PEST ++ user manual version 5.2.18. Github. Available online at: https://github.com/usgs/pestpp/blob/master/documentation/pestpp_users_manual.md (Accessed April 23, 2025).

Google Scholar

Puckett, L. J., and Cowdery, T. K. (2002). Transport and fate of nitrate in a glacial outwash aquifer in relation to ground water age, land use practices, and redox processes. J. Environ. Qual. 31 (3), 782–796. doi:10.2134/jeq2002.7820

PubMed Abstract | CrossRef Full Text | Google Scholar

Reed, E. M., Wang, D., and Duranceau, S. J. (2018). Evaluating nitrate management in the volusia blue springshed. J. Environ. Eng. 144 (3), 05018001. doi:10.1061/(ASCE)EE.1943-7870.0001324

CrossRef Full Text | Google Scholar

Schullehner, J., Hansen, B., Thygesen, M., Pedersen, C. B., and Sigsgaard, T. (2018). Nitrate in drinking water and colorectal cancer risk: a nationwide population-based cohort study. Int. J. Cancer 143 (1), 73–79. doi:10.1002/ijc.31306

PubMed Abstract | CrossRef Full Text | Google Scholar

Shakeri, R., Nassery, H. R., and Ebadi, T. (2023). Numerical modeling of groundwater flow and nitrate transport using MODFLOW and MT3DMS in the karaj alluvial aquifer, Iran. Environ. Monit. Assess. 195, 242. doi:10.1007/s10661-022-10881-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Shamrukh, M., Corapcioglu, M. Y., and Hassona, F. A. (2001). Modeling the effect of chemical fertilizers on ground water quality in the nile valley aquifer, Egypt. Groundwater 39 (1), 59–67. doi:10.1111/j.1745-6584.2001.tb00351.x

CrossRef Full Text | Google Scholar

Sieczka, A., Bujakowski, F., Falkowski, T., and Koda, E. (2018). Morphogenesis of a floodplain as a criterion for assessing the susceptibility to water pollution in an agriculturally rich valley of a lowland river. Water 10 (4), 399. doi:10.3390/w10040399

CrossRef Full Text | Google Scholar

Smith, E. A., and Westenbroek, S. M. (2015). Potential groundwater recharge for the state of Minnesota using the soil-water-balance model, 1996-2010. U.S. Geological Survey Scientific Investigations Report 2015-5028. doi:10.3133/sir20155038

CrossRef Full Text | Google Scholar

Stark, J. R., Armstrong, D. S., and Zwilling, D. R. (1994). “Stream-aquifer interactions in the straight river area,” in Becker and hubbard counties. Minnesota: US Geological Survey Water-Resources Investigations Report, 94–4009. doi:10.3133/wri944009

CrossRef Full Text | Google Scholar

Stroom, K. (2024). “Straight river (hubbard county, Minnesota) nutrient study,” in Minnesota pollution control agency. Available online at: https://www.pca.state.mn.us/sites/default/files/wq-ws5-07010106c.pdf (Accessed October 23, 2024).

Google Scholar

Tajtáková, M., Semanová, Z., Tomková, Z., Szökeová, E., Majoroš, J., Rádiková, Ž., et al. (2006). Increased thyroid volume and frequency of thyroid disorders signs in schoolchildren from nitrate polluted area. Chemosphere 62 (4), 559–564. doi:10.1016/j.chemosphere.2005.06.030

PubMed Abstract | CrossRef Full Text | Google Scholar

United States Department of Agriculture, National Agriculture Statistics Service (USDA NASS) (2025). Crop data layer. Available online at: https://nassgeodata.gmu.edu/CropScape/(Accessed October 5, 2025).

Google Scholar

United States Department of Agriculture, National Agriculture Statistics Service (USDA NASS) (2024). Potatoes 2023 summary. Available online at: https://downloads.usda.library.cornell.edu/usda-esmis/files/fx719m44h/6q184c88z/v118t7605/pots0924.pdf (Accessed April 24, 2025).

Google Scholar

United States Department of Agriculture, Risk Management Agency (USDA RMA) (2024). PM-24-002: 2023 crop year (CY) dry bean and dry pea revenue endorsement harvest prices. Available online at: https://old.rma.usda.gov/en/Policy-and-Procedure/Bulletins-and-Memos/2024/PM-24-002 (Accessed October 5, 2025).

Google Scholar

United States Environmental Protection Agency (US EPA) (2007). Nitrates and nitrites: TEACH chemical summary. U.S. EPA toxicity and exposure assessment for children’s health (TEACH) database. Available online at: https://archive.epa.gov/region5/teach/web/pdf/NO3-Ns_summary.pdf (Accessed October 23, 2024).

Google Scholar

Ward, M. H., Kilfoy, B. A., Weyer, P. J., Anderson, K. E., Folsom, A. R., and Cerhan, J. R. (2010). Nitrate intake and the risk of thyroid cancer and thyroid disease. Epidemiology 21 (3), 389–395. doi:10.1097/EDE.0b013e3181d6201d

PubMed Abstract | CrossRef Full Text | Google Scholar

Ward, M. H., Jones, R. R., Brender, J. D., De Kok, T. M., Weyer, P. J., Nolan, B. T., et al. (2018). Drinking water nitrate and human health: an updated review. Int. J. Environ. Res. Public Health 15 (7), 1557. doi:10.3390/ijerph15071557

PubMed Abstract | CrossRef Full Text | Google Scholar

Wayment, J. (2021). Nitrate leaching mitigation with kura clover and rye covers for corn and soybean in irrigated sands. St. Paul (MN): University of Minnesota.

Google Scholar

White, J. T. (2018). A model-independent iterative ensemble smoother for efficient history-matching and uncertainty quantification in very high dimensions. Environ. Model. and Softw. 109, 191–201. doi:10.1016/j.envsoft.2018.06.009

CrossRef Full Text | Google Scholar

White, J. T., Knowling, M. J., Fienen, M. N., Feinstein, D. T., McDonald, G. W., and Moore, C. R. (2020). A non-intrusive approach for efficient stochastic emulation and optimization of model-based nitrate-loading management decision support. Environ. Model. and Softw. 126, 104657. doi:10.1016/j.envsoft.2020.104657

CrossRef Full Text | Google Scholar

White, J. T., Knowling, M. J., Fienen, M. N., Siade, A., Rea, O., and Martinez, G. (2022). A model-independent tool for evolutionary constrained multi-objective optimization under uncertainty. Environ. Model. Softw. 149, 105316. doi:10.1016/j.envsoft.2022.105316

CrossRef Full Text | Google Scholar

Wilson, M. L., Moncrif, J. F., and Rosen, C. J. (2008a). Kidney bean (Phaseolus vulgaris L.) production on an irrigated, coarse-textured soil in response to polymer coated urea and tillage: II. Plant N accumulation, nitrate leaching and residual inorganic soil N. J. Environ. Monit. Restor. 5 (1), 77–91. doi:10.4029/2008jemrest5no18

CrossRef Full Text | Google Scholar

Wilson, M. L., Moncrif, J. F., and Rosen, C. J. (2008b). Kidney bean (Phaseolus vulgaris L.) production on an irrigated, coarse-textured soil in response to polymer coated urea and tillage: I. Grain yields, disease severity, and a simple economic analysis. J. Environ. Monit. Restor. 5 (1), 41–56. doi:10.4029/2008jemrest5no15

CrossRef Full Text | Google Scholar

Yang, L., Zheng, C., Andrews, C. B., and Wang, C. (2021). Applying a regional transport modeling framework to manage nitrate contamination of groundwater. Groundwater 59 (2), 292–307. doi:10.1111/gwat.13047

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, C., and Bennett, G. D. (2002). Applied contaminant transport modeling. 2nd ed. New York: Wiley-Interscience.

Google Scholar

Zheng, C., and Wang, P. P. (1999). MT3DMS: a modular three-dimensional multispecies transport model for simulation of advection, dispersion, and chemical reactions of contaminants in groundwater systems; documentation and user's guide. Washington, DC: U.S. Army Corps of Engineers Engineer Research and Development Center Contract.

Google Scholar