Eutrophic conditions in estuaries in the province of Prince Edward Island (PEI), Canada have been documented since 2002 with the implementation of the PEI Estuaries Survey (Bugden et al., 2014), although they have been anecdotally noted for decades prior. PEI is intensively cultivated for potatoes, accounting for 25% of Canadian potato exports (Government of Prince Edward Island, 2020). Potatoes require high fertilization rates but have poor nitrogen (N) uptake efficiency (Delgado et al., 2001), which causes significant leaching and runoff of N into surface and groundwater (De Notaris et al., 2018; Jiang et al., 2011; Zebarth et al., 2009). The island is characterized by well-drained soils and shallow aquifers that are highly susceptible to N contamination in the form of NO₃-N (Puckett et al., 2011). Rainfall, snowmelt, and associated recharge elevate water tables in PEI and drive discharge of N-rich groundwater to streams and, eventually, PEI coastal waters (Pavlovskii et al., 2023). Under current agricultural practices, groundwater NO₃-N concentrations in PEI are predicted to continue to increase as NO₃-N penetrates deeper into bedrock aquifers (Jiang and Somers, 2009; Zhang and Hiscock, 2011). This could deleteriously impact potable water quality in this province, which is fully dependent on groundwater for drinking water supply (Council of Canadian Academies, 2009), as well as the health of coastal ecosystems that receive contaminated groundwater.

Amidst growing public concern surrounding the health of PEI’s environment, the province has implemented some mandatory agricultural BMPs since the early 2000s. For example, 15 m buffer strips are required around any watercourse (Government of PEI, 2016), and the Agricultural Crop Rotation Act (Government of PEI, 2002) states that potatoes can only be grown on the same parcel of land once every 3 years. The Province of PEI also provides funding through the Agriculture Stewardship Program for farmers wishing to implement BMPs under the categories of manure and livestock management, agroforestry, water and supply management and integrated pest management (Sustainable CAP, 2023). Examples of BMPs currently employed in PEI are cover cropping, incorporating N-fixing plants into crop rotations, improving liquid manure application, implementing edge of field treatment systems such as buffer strips, constructed wetlands and vegetative swales. Research into the effects of alternative crop rotations to the traditional Potato-Barley-Clover (PBC) rotation on PEI has revealed that replacing red clover with soybeans, a legume, adds an additional cash crop, while reducing NO₃-N leaching and increasing potato tuber yield (Azimi et al., 2022; Liang et al., 2019).

Basin Head lagoon, located on the northeastern coast of PEI (Figure 1A), has been experiencing declining ecosystem health (Connolly, 2002). Accordingly, the lagoon was designated as a federal Marine Protected Area in 2005 with the goal of protecting the endemic giant Irish moss, a unique morphotype of Irish moss (Chondrus crispus), from declining ecosystem health and water quality (DFO, 2009). More prevalent and sustained hypoxic events have been noted over recent years, particularly in the shallower and more poorly mixed northeast arm of the lagoon throughout the spring and summer months (DFO, 2021). Hypoxic or anoxic conditions can result in irreversible changes to aquatic community structure (Coffin et al., 2018) and can threaten the health of the Irish moss in the Basin Head Marine Protected Area.

The Basin Head watershed is 14.6 km² (Figure 1B) with agriculture (41%) and forest (31%) comprising most of the land cover (KarisAllen et al., 2022). This study focuses on the four major tributaries that discharge into the lagoon with a collective drainage area of 7.2 km² (Figure 1B), and does not include direct runoff to the lagoon. The Charlottetown soil series is the most abundant throughout the entire province (MacDougall et al., 1988), and is the dominant soil type in the Basin Head watershed (95%). It is classified as a well-drained sandy loam located on gentle slopes (4%–9%) and is extremely well suited for agriculture (MacDougall et al., 1988). The porosity of the Charlottetown soil series varies from 30% in the till layer below 0.5 m (Francis, 1989; Heath, 1983; MacDougall et al., 1988) to approximately 50% in the upper, macroporous zone (Carter, 1987; MacDougall et al., 1988). Depth to bedrock, based on historical well data in the Kings County region, has been observed to range from 0 to 6.2 m (Government of PEI, 2019; Joostema, 2015) and was measured to be 4.6 m at a shallow monitoring well in the watershed installed near the lagoon (KarisAllen et al., 2022) (Figure 1B). The bedrock aquifer that underlies PEI is primarily comprised of sandstone (80%–85%) and interbedded mudstone which is fractured, creating preferential flow paths (Francis, 1989). The fractured sandstone aquifer functions as a dual-porosity system, which presents challenges with characterization of aquifer effective porosity, specific yield, and hydraulic conductivity. However, equivalent porous medium approaches are commonly applied when modeling PEI groundwater systems (e.g., Stanic et al., 2024).

The climate of PEI is characterized as temperate and humid, with a mean annual precipitation of 1,173 mm at the Charlottetown weather station (1971–2000) and mean annual air temperature of 5.3°C (Paradis et al., 2016). Most precipitation (75%) falls as rain (Paradis et al., 2016), and precipitation is relatively uniformly distributed throughout the year, averaging between 80 and 120 mm per month (ECCC, 2023). Baseflow typically contributes 60%–70% of streamflow annually on PEI, and up to 100% in the summer months (Benson et al., 2007; DELJ, 2013); hence groundwater is critical to quantifying PEI water balances and nutrient transport. The fractured sandstone aquifer in the Basin Head watershed provides baseflow to the four tributaries that are the focus of this study (Figure 1) and provides flow for the 30+ intertidal springs that discharge directly into the lagoon (KarisAllen and Kurylyk, 2021; Oliver et al., 2024). The lagoon is tidally pumped and experiences mixing from waves, surges, and current during coastal storms (Bonnington et al., 2023).

The SWAT+ model (Arnold et al., 2012) requires climate variable inputs as well as GIS-based inputs for topography (i.e., Digital Elevation Model, DEM), land use and soils. A summary of the spatial data used to construct the SWAT+ model is presented in Table 1. Full maps of DEM, soil, and land use are provided in Supplementary Figures S1–S3. Figure 2 displays all soil types within the watershed boundaries; however, for simplicity in the model, the Charlottetown series was chosen to represent the entirety of the watershed due to its relative abundance (95%). Precipitation and maximum and minimum daily air temperature records were obtained from a weather station (Onset HOBO) installed at Basin Head (Figure 1B). This weather station did not provide accurate precipitation in winter months (December to March) as it did not possess a heated rain gauge; therefore, the record was supplemented with values from the nearest Environment and Climate Change Canada (ECCC) weather stations at East Point (ECCC ID 7177) and St. Peters (ECCC ID 7177), which are 12 km and 36 km from Basin Head. The East Point station was used to develop a local weather generator within the model (Neitsch et al., 2011) to simulate other variables such solar radiation, wind speed and relative humidity. The weather generator uses statistical monthly measures for air temperature, precipitation, solar radiation, dew point, and more to develop representative daily climate data for a subbasin.

SWAT+ models three major forms of nitrogen in the soil profile and shallow aquifer: (1) organic nitrogen associated with humus, (2) mineral forms of nitrogen held by soil colloids, and (3) mineral forms of nitrogen in solution. The movement of nitrogen between two ‘pools’ (i.e., valence states) of inorganic nitrogen, ammonium (NH₄ ⁺) and nitrate (NO₃ ⁻), and three pools of organic nitrogen is simulated. Nitrate transport in SWAT+ is simulated only in the mobile water layer, which includes water transported by surface runoff, lateral flow (i.e., groundwater flow within the soil profile) or percolation (i.e., water moving past the lowest soil profile layer and into the shallow aquifer). The mobile water layer nitrate concentration, c o n c N O 3 , m o b i l e , is multiplied by the flow from each pathway (surface runoff, lateral flow and percolation) to determine nitrate load from each source in the soil layer. Important parameters in these partitioning equations are θ e , the fraction of porosity where anions are excluded, and S A T l y , the saturated water content of the soil layer (mm H₂O), which determine nitrate load from each source in the soil layer, and δ g w , the delay time associated with overlying geologic formations (days), which determines nitrate load in the shallow aquifer.

Another important input to the SWAT+ model is agricultural practices. Fertilizer amount, type, and application timing, as well as crop rotation timing influences the forcing of chemical leaching into soils and surface water runoff, and hence NO₃-N concentrations in receiving water bodies. Several groups with local knowledge of agricultural practices in Kings County, PEI, where Basin Head is located (i.e., Souris Wildlife, PEI Potato Board, and Fisheries and Oceans Canada) were contacted to develop representative agricultural practices for model inputs. Cropland in the watershed is assumed to adhere to local industry standard crop rotations of potato, a grain (e.g., barley), and a forage (e.g., clover) on a 3-year rotation (herein PBC) as is mandated by the province (Azimi et al., 2022). Although in recent years clover has been replaced with a mix of legumes including alfalfa, sudangrass, pearl millet, and ryegrass, red clover was simulated to maintain consistency over the period of interest (1990–2100). Standard fertilizer application to potato crops in PEI is 155 kg N ha⁻¹ (PEI Analytical Laboratories, Department of Agriculture and Fisheries, PEI). A 15-15-15 Nitrogen-Phosphorous-Potassium fertilizer was chosen based on availability in the SWAT+ interface and given the similarity to the 17-17-7 fertilizer employed by Liang et al. (2020) who studied N dynamics in PEI under different crop rotations. The application method of broadcast, which refers to fertilizer being spread uniformly over the entire field, was chosen based on local practices in the Basin Head watershed, Supplementary Table S1 provides a summary of the agricultural operations incorporated in SWAT+.

Discrete water samples were collected at the four main (i.e., highest flows) tributaries on 15 occasions covering diverse flow conditions from 11 November 2021, to 28 November 2023 (Figure 1B). Samples were analyzed for NO₃-N using the cadmium reduction colorimetric method-APHA Method 4500-E (APHA, 2022). These grab samples were also analyzed for total suspended solids (TSS) using Standard Method 2450D (APHA, 2022). Continuous stream stage measurements were recorded using Onset HOBO loggers at 15-min intervals. Stage was barometrically compensated with a HOBO air pressure logger (Figure 1B). Instantaneous stream discharge measurements were obtained using the velocity-area method (Turnipseed and Sauer, 2010) during the period of 26 June 2019, to 14 October 2023. Stage-discharge relationships were developed for each tributary (Supplementary Figure S4) and used to convert stage to flow. Monthly streamflow (Supplementary Figure S5) was used to calibrate the SWAT+ model hydrology, which is a standard interval for such SWAT+ applications (Ricci et al., 2023; Van Liew et al., 2012), while water quality measurements on a daily basis (i.e., sample concentration) were used to calibrate SWAT+ parameters for NO₃-N and TSS.

In groundwater-dominated watersheds, generally the most sensitive streamflow model parameters relate to groundwater and soil moisture. A literature review for SWAT+ modeling found the SCS curve number (cn2), baseflow recession constant (alpha), fraction of root zone percolation that reaches deep aquifer (perco), and soil evaporation compensation factor (esco) to be most sensitive (Ahmad et al., 2011; Liang et al., 2020; Spruill et al., 2000). These parameters, along with others identified in a sensitivity analysis conducted in SWAT+ Toolbox (Supplementary Table S2), were chosen to calibrate streamflow (James, 2022). Evaluation statistics used to assess model performance included the Nash-Sutcliffe efficiency (NSE), RMSE-observation standard deviation ratio (RSR), and percent bias (PBIAS) (Moriasi et al., 2007). Moriasi et al. (2007) states that for streamflow, an RSR less than 0.7 is satisfactory, an NSE over 0.4 is satisfactory, and PBIAS less than 10 is very good. The corresponding equations are presented below, where Y_i ^obs is the observed streamflow and Y_i ^sim is the simulated streamflow generated by SWAT+. NSE = 1 − ∑ i = 1 n Y i obs −   Y i sim 2 ∑ i = 1 n Y i obs −   Y i mean 2 (1) RSR = ∑ i = 1 n Y i obs −   Y i sim 2 ∑ i = 1 n Y i obs −   Y i mean 2 (2) PBIAS = ∑ i = 1 n Y i obs −   Y i sim × 100 ∑ i = 1 n Y i obs (3)

Grab sample concentrations (Supplementary Table S3) were used to calibrate NO₃-N outputs from SWAT+. Daily NO₃-N mass (kg/day) output from each tributary from SWAT+ were converted to a concentration (mg/L) by dividing the mass by the daily total flow. Modeled concentrations were compared to 11 grab sample concentrations from 11 November 2021, to 5 July 2023, and Root Mean Squared Error (RMSE) was calculated for the modeled vs. measured concentrations: RMSE = ∑ i = 1 n y i ^ − y i 2 n (4) where y n ^ are predicted values, y n are observed values, and n is the number of observations. RMSE was used as the statistical performance metric for NO₃-N and TSS calibration as continuous measurements were not available for these parameters. A SWAT study conducted on the relatively nearby Wilmont River watershed in southwestern PEI (Liang et al., 2020) was used as a reference for the calibration in the present study due to the similarity in study region and hydrologic conditions. Supplementary Table S4 shows the parameters used in the sensitivity analysis and calibration in Liang et al. (2020), which were chosen as a starting point for the calibration of NO₃-N at the study site.

Once the model is calibrated satisfactorily, projected results using climate change data can be analyzed to determine patterns in seasonality of flow and NO₃-N loading, and trends to baseflow index (BFI). BFI is calculated by dividing the total volume of baseflow by total streamflow for a period. The SWAT+ outputs of lateral flow (flo) and shallow groundwater flow (latQ) sum to determine baseflow, and total streamflow is the summation of water yield (wateryld) and lateral flow (flo) into the channel. The ‘half-flow date’, T50, is a measure of streamflow timing and is defined as the date by which 50% of streamflow for a given water-year (i.e., October 1 to September 30) has passed (Court, 1962). This measure has been used to characterize flow regimes in the neighboring province of Nova Scotia (Johnston et al., 2022), and was calculated in this study for the three model periods.

It should be noted that SWAT studies in other small (<10 km²) watersheds have been noted to yield unsatisfactory results (Spruill et al., 2000), largely due to characteristically low times of concentration and a propensity toward ‘flashiness’ in peak runoff (Uzeika et al., 2012). Chu et al. (2004) noted that characterizing baseflow correctly can be more important in small watersheds. Bailey et al. (2020) outlines limitations to the current groundwater algorithms in SWAT+ and proposes an amendment to the standard SWAT+ module to improve model accuracy, particularly in baseflow-dominated watersheds where model performance can be poor (Bosch et al., 2010; Luo et al., 2012). Despite the challenges in applying SWAT+ in small, baseflow-dominated systems, this model is widely used and still provides useful insight into the hydrology and transport dynamics in these environments.

Changes to climate forcing will drive changes in hydrologic processes and agricultural practices (Gordon et al., 2008), which can be represented in SWAT+. An analysis by the Pacific Climate Impacts Consortium [of 26 downscaled Coupled Model Intercomparison Project (CMIP) 6 Global Climate Models (GCMs)] (Eyring et al., 2016) and three Shared Socioeconomic Pathways (SSPs) (Riahi et al., 2017) was completed to select models projecting low, moderate, and high changes to precipitation and temperature. Climate data for these models were downloaded from the NASA Climate Data Services (NASA, 2023), and, for the sake of downscaling consistency, only models downscaled using Bias Correction/Constructed Analogues with Quantile Mapping Reordering (BCCAQv2) (Maurer et al., 2010; Cannon, 2015) were considered. The scenarios from the IPSL-CM6A-LR and EC-Earth3-veg GCMs were removed from consideration because they are deemed ‘hot’ models, projecting much greater warming than is generally expected (Rahimpour Asenjan et al., 2023). UKESM1-0-LL was removed for simplicity since it only represents a 360-day year, and BCC-CSM2-MR was missing from the NASA dataset. The remaining six models were compared based on their average annual precipitation for each modeled time period (1990-2020, 2040-2070 and 2070–2100) (Supplementary Table S5), and three representative models (INM-CM5-0, FGOALS-g3, and MRI-ESM2-0) were chosen for low, moderate, and high climate change scenarios (Figure 2). Two SSPs per climate model, SSP2-4.5 (RCP 4.5) and SSP5-8.5 (RCP 8.5), were chosen to represent mid-range and high emissions scenarios, yielding a total of six climate change scenarios with a wide range of possible climate outcomes. The climate data representing each scenario was entered in SWAT+, and the calibrated model was run from 1990 to 2100, with outputs from representative periods of historical (1990–2020), mid-century (2040–2070) and end of century (2070–2100) used for further analysis. The observation period GCM results for each climate scenario were used for the 1990-2020 comparison period rather than historic data in the watershed because the climate station in the watershed was only operational since 2018, and there is no nearby long-term climate station or one with very similar measurements to the Basin Head climate station (Figure 1B) during its operation. Therefore, the climate change analysis presented herein should be considered a sensitivity analysis to plausible future projections of climate change when compared to historic-period output for the same downscaled GCMs. A one-way ANOVA was conducted in Minitab (2021) to generate monthly pairwise comparisons of the NO₃-N loading between each period (N = 30) to determine if results were significantly different (p < 0.0.5) among scenarios.

It is important to consider societal changes that may arise due to climate change. Climate change will result in extended growing seasons in Atlantic Canada, defined by the number of days per year where the air temperature is over 5°C, with a projected increase in growing season length of 11% by 2050 and 18% by 2080 (Richards and Daigle, 2011) using climate data from the nearest station to Basin Head (East Baltic, 46.43N 62.17W). In SWAT+, these percentages were rounded to a 10% increase in growing season length for the mid-century period (2040–2070) and a 20% increase for the end of century period (2070–2100). To account for a longer growing season and potential loss in plant productivity, an additional percentage of fertilizer applied to crops was also incorporated. Arbitrarily, an increase in fertilizer application of 10% and 20% were applied for mid-century and end of century periods, respectively, to align with the projected increase in growing season length (Richards and Daigle, 2011). This future agricultural scenario was only run for the MRI-ESM2-0 SSP5-8.5 GCM to consider changing agricultural practices for a more extreme climate scenario. All model runs are listed in Supplementary Table S6 with corresponding weather inputs.

BMPs are implemented widely across PEI and include crop rotations, buffer strips, cover crops, farmable berms, and terraces. Three common BMPs were modeled in SWAT+ from 1990 to 2100 under the climate change scenario with the highest cumulative precipitation (MRI-ESM2-0 SSP5-8.5) to understand how NO₃-N loading could be mitigated in a ‘worst-case’ climate scenario. The BMP scenarios chosen to be simulated in SWAT+ were (1) buffer strips (60 m, 30 m), (2) delayed tillage of forage crop from fall to spring, and (3) alternative crop rotation of Potato-Soybean-Barley (PSB) from the default PBC rotation used for other SWAT+ runs.

Buffer strips refer to narrow strips of dense grass surrounding waterways on agricultural land (Helmers et al., 2006). They are known to slow down runoff velocity to control field erosion and trap particulate pollutants and have been shown to facilitate NO₃-N removal from shallow groundwater (Simpkins et al., 2002). A spatial analysis revealed that most agricultural fields in the watershed have a grassed buffer of 30 m between the field and waterways, which is double the provincially mandated width of 15 m. A scenario incorporating a 60 m buffer was simulated, which would be considered an extreme measure. The removal of NO₃-N in vegetative filter strips in SWAT+ is predicted based on the ratio between field area and filter strip area. It should be noted that the decrease in cropland associated with an increase in buffer size was not included in this simulation to isolate the effect of the buffer, thus, the total effect of an increased buffer width would be underestimated in these results.

Typically, on PEI, a forage crop (i.e., clover, timothy, rye) is grown in the year prior to potato planting. By delaying tillage of forage from fall to spring, just before planting, NO₃-N leaching in the winter months is reduced, and there is residual NO₃-N left in the soil in spring, reducing the need for fertilizer (Stuart, 2017). This effect was modeled in SWAT+ by editing the agricultural operations such that clover was plowed May 15th, prior to planting of potatoes, instead of on September 15th and by reducing fertilizer application by 25% (Stuart, 2017).

PBC is a representative crop rotation for PEI; however, recent studies have suggested that incorporating legumes may improve soil health, increase potato yields and reduce NO₃-N surplus leading to leaching due to the N fixing characteristics of legumes (Liang et al., 2019; Whittaker et al., 2023). A PSB crop rotation was modeled in SWAT+ (Supplementary Table S7) to test this hypothesis under the three periods, for which fertilizer application was reduced by 50% to account for the N fixing characteristics of soybeans. This reduction amount was chosen arbitrarily based on the results of Liang et al. (2019) that showed N efficiency (i.e., how well plants utilize available N) of PSB crop rotations were 1.6 times that of PBC rotations.

73 water samples from the Basin Head tributaries were analyzed for NO₃-N, of which 24 (33%) exceeded the Canadian Council of Ministers of the Environment (CCME) guideline for the protection of aquatic life in freshwater (3 mg L⁻¹, CCME, 2012). NO₃-N concentrations range from 0 to 6.2 mg/L, with tributary 1 having the highest average concentration (4.2 mg/L) compared to the other tributaries (1.3–1.9 mg/L). TSS concentrations ranged from 0 to 44.0 mg/L, with tributary 2 and 3 being higher on average (11.3–12.3 mg/L) than tributary 1 and 2 (6.1–7.2 mg/L). All grab sample measurements are summarized in Supplementary Table S3.

All sensitive model parameters for hydrology (Supplementary Table S8) were manually calibrated to the measured streamflow (Supplementary Figure S5) until the goodness of fit was satisfactory from April 1st to 1 September 2022, to represent hydrology in the growing season (Table 2). It should be noted that model results from winter months are not calibrated, thus one must take caution when drawing conclusions during these months. Also, snowmelt_min, snowmelt_max, and snowmelt_tmp are climate-related parameters, and would most likely be spatially uniform for a small watershed; they were manually changed to improve streamflow calibration. The sensitivity analysis revealed that simulated NO₃-N concentrations were sensitive to only six out of eight parameters: hlife_n, the half-life of nitrate in the shallow aquifer, nperco, the nitrate percolation coefficient, sdnco, the denitrification threshold water content, cdn, denitrification exponential rate coefficient, cmn, the rate factor for humus mineralization of active organic nutrients and n_updis, the nitrogen uptake distribution parameter. NO₃-N and TSS were calibrated from 11 November 2021, to 5 July 2023, and from 27 July 2022, to 5 July 2023, respectively. The graphical representation of NO₃-N calibration and RMSE (Equation 4) values between measured to simulated concentration are found in Supplementary Figure S6. Six parameters were calibrated for NO₃-N: hlife_n (600 for tributary 1, 200 for tributaries 2–4), nperco (0.15 for tributaries 1-3, 0.05 for tributary 4), sdnco (0.66), cdn (0.07), cmn (0.0029) and n_updis (66.7). RMSE (Equation 4) between measured NO₃-N load (i.e., concentration multiplied by flow) and simulated load were calculated to be 2.7, 7.6, 4.5 and 1.0 kg/day for tributaries 1 through 4.

The sediment calibration process in this study closely followed a previous SWAT study (Sinclair, 2014) by choosing sensitive starting parameters for calibration with a focus on in-channel processes. SWAT+ uses the Modified Universal Soil Loss Equation (MUSLE) to generate the landscape erosion contribution to channel TSS values. MULSE is a function of several factors, which were entered into SWAT+ as default values, except for the soil erodibility factor, which was calculated to be 0.49 based on soil information. High sediment loads in several PEI rivers have been documented, mainly citing sediment deposition resulting from intensive agricultural practices (Sirabahenda et al., 2020). Supplementary Figure S7 presents a graphical representation of TSS calibration and RMSE values between measured and simulated concentrations. Four parameters were calibrated for TSS: spexp (1), spcon (0.0001), cov fact (0.02 for tributaries 2 and 3, 0.005 for tributaries 1 and 4), and erod_fact (0.25 for tributaries 1 and 4, 0.75 for tributary 2 and 0.5 for tributary 3). RMSE between measured TSS load (i.e., concentration multiplied by flow) and simulated load were calculated to be 6.7, 34.7, 36.5, and 1.1 kg/day for tributaries 1 through 4. Final calibrated RSME values for NO₃-N, TSS, and streamflow for each tributary are found in Table 2. NO₃-N was successfully calibrated; however, RMSE values for TSS are significantly higher than 1–2 mg/L, indicating that it is only preliminarily calibrated for sediment transport.

Modelled conditions using MRI-ESM2-0 SSP5-8.5 from the historical period (1990–2020) generally aligned with expected water balance for PEI where an annual average of 1,100 mm of precipitation is partitioned into 440 mm evapotranspiration, 260 mm surface water runoff and 400 mm is groundwater flow (PEI Department of Fisheries and Environment, 1996). During this period SWAT+ simulates the average annual NO₃-N transported via surface runoff as 1 kg/ha, while average annual NO₃-N transported through groundwater is 6 kg/ha, which aligns with other PEI studies such as Grizard. (2013) that simulated an average NO₃-N load of 12 kg-N/ha watershed/yr for the period 1996-2012. Average annual sediment load via overland erosion was modelled to be 0.96 t/ha/yr, which aligns with a nearby previous study in Souris River, PEI that found overland erosion to occur at a rate of 1.5–1.8 t/ha/yr, attributing this high sediment load to row crop agriculture which is common across PEI (Alberto et al., 2016). SWAT+ models deposition of sediment in all tributaries in the watershed, meaning that all TSS at the channel outlets are originated from overland erosion.

Changes to moisture conditions and seasonality of precipitation events due to climate change could result in changes to the ratio of baseflow and surface runoff contributions in streamflow. The monthly averaged BFI was analyzed for each period in SWAT+, revealing shifts in the partitioning of streamflow over time (Figure 4B). Average BFI in the months of January to May increased from 0.39 in the historical (1990–2020) period to 0.74 and 0.78 in the mid-century (2020–2070) and end-of-century (2070–2100) periods. The largest BFI change occurred in February where monthly average BFI increased from 0.12 in the historical period to 0.63 and 0.68 in the two future periods (>five-fold increase). Average annual simulated BFI increased 22% from historical (1990–2020) to mid-century (2040–2070) periods, and only increased 1% between mid-century and end of century (2070–2100) periods. These baseflow changes could be attributed to increased precipitation as rainfall, more mid-winter thaws, and an earlier spring melt period, inducing recharge and subsequent flow through the shallow aquifer during the winter season.

The SWAT+ model projections also show suppressed summer flow, particularly in April and May, and increased winter flow between historical (1990–2020) and both future periods (2040-2070 and 2070–2100) and illustrate that projected changes in nitrate loading and streamflow follow similar seasonal patterns (Figures 4A, C). During the historical period in the months of January and February, the standard deviation in flow is much lower than the end-of-century period (average decrease of 74%). However, in the months of April to May, the standard deviations of the flow values in the end-of-century period are much lower than in the historical period (−72% on average). This indicates that in the winter months (January and February), there is more variance in flow in the end-of-century period, while in the spring months (April and May) there is more variance during the historical period. Standard deviations of flow between historical and both future periods were within 50% of each other from June to November and were highest from December to May. Increased flow variability in winter and early spring could reflect uncertainty in the effects of climate change on snowmelt processes. The historical period T50 (March 19), a measure of streamflow timing, was shifted earlier by 73 and 78 days respectively for the mid-century (January 7) and end of century (January 1) periods. Some studies (Mukundan et al., 2020; Shrestha et al., 2012) have proposed that earlier snowpack melt and an increase to annual rainfall due to climate change, both of which were modeled in this study, will increase the magnitude and shift the timing of nutrient loading. Average annual streamflow is simulated to increase 5% between historical (1990–2020) and mid-century (2040–2070) periods, and increase 4% between mid-century and end of century (2070–2100) periods. This could be attributed to increased precipitation expected in a future climate, or increased snowmelt due to higher temperatures. Higher streamflows with on average higher baseflow influence in future years (2040–2100) could result in a higher flux of NO₃-N to coastal waters.

Based on the SWAT+ simulations of future climate change impacts and the corresponding ANOVA analysis, only the scenario with the highest cumulative precipitation (MRI-ESM2-0 SSP5-8.5) resulted in statistically significant (p < 0.05) changes in NO₃-N loading to the lagoon, and only when comparing historical to end-of-century (1990-2020 to 2070–2100) time periods (Figure 3). Under this scenario, the median annual NO₃-N load increased from 5,752 to 7,239 kg (+26%). Among all climate scenarios considered, this scenario produced both the highest mean annual precipitation and the highest mean annual air temperature at the end of century (Supplementary Figures S8, S9).

The GCM scenarios with the lowest cumulative precipitation (INM-CM5-0 SSP2-4.5 and INM-CM5-0 SSP5-8.5) saw small changes in median NO₃-N load of less than 7% between consecutive model periods. The mid-range cumulative precipitation GCM scenario, FGOALS-g3 SSP2-4.5, yielded interesting results, as annual median NO₃-N load increased by 13% between historical (1990–2020) and mid-century (2040–2070) periods but decreased by 2% from mid-century (2040–2070) to end of century (2070–2100). The FGOALS-g3 SSP5-8.5 GCM yielded a 5% and 10% increase in annual median NO₃-N load between the same periods. In general, the simulations with higher cumulative precipitation resulted in larger increases in annual NO₃-N load (Figure 2). Between each subsequent modeled period, the MRI-ESM2-0 SSP2-4.5 GCM resulted in a 9% and 6% increase to median annual NO₃-N load, while MRI-ESM2-0 SSP5-8.5 resulted in an increase of 13% and 11%, respectively.

These NO₃-N load trends generally correspond to trends in cumulative precipitation (Figure 2) and suggest that the projected increased NO₃-N loading could be attributed to increased water yield and/or enhanced biogeochemical cycling of N (e.g., greater mineralization of organic N and nitrification of ammonia-N) (Pesce et al., 2018; Rabalais et al., 2009). The similarity in pattern of seasonal streamflow and NO₃-N (Figures 4A, C) indicates that increased water yield, rather than elevated concentration, is predominantly responsible for the projected increase in NO₃-N loading.

The effects of NO₃-N leaching through surface and groundwater pathways under future climate change scenarios is presently not well understood. However, it is known that climate change impacts groundwater recharge and discharge (Kurylyk et al., 2014), sea levels, soil conditions (i.e., soil moisture, organic carbon, alkalinity) and agricultural productivity (Stuart et al., 2011), all of which could impact nutrient loading and transport dynamics in coastal watersheds. Some studies investigating the effects of climate change on nutrient loading highlight the uncertainty of climate change models on a regional scale (Bürger et al., 2013; El-Khoury et al., 2015; Jeppesen et al., 2011). These have suggested that the use of different downscaling methods for GCMs can result in conflicting trends in directions and magnitudes of nutrient loading, which is likely partly related to the uncertainty in future trend directions and magnitudes for precipitation and groundwater recharge (Kurylyk and MacQuarrie, 2013). A few studies (e.g., Luo et al., 2020; Qi et al., 2009; Tong and Liu, 2006) have suggested that changes to hydrology have more influence on nutrient loading dynamics than changes in land use, while a local PEI study (De Jong et al., 2008) suggests land use is the more important factor. Pulido-Velazquez et al. (2015) noted the importance of analyzing the impacts of land use and hydrologic stressors together when considering nutrient loading patterns. The use of SWAT+ in this study addresses many of these concerns by incorporating future climate data and current land use and agricultural practices, with capability to simulate a wide range of future outcomes.

Seasonal trends in NO₃-N loading to the lagoon were also analyzed under the highest cumulative precipitation climate scenario, MRI-ESM2-0 SSP5-8.5 (Figure 4C). Simulations during the historical period (1990–2020) indicate that NO₃-N loading to the lagoon is highest during the months of May and November. These peaks likely correspond to NO₃-N availability due to recent fertilizer application in spring or harvest in fall and corresponding large precipitation events and snowmelt that typically occur in spring and fall in PEI. For the future climate scenarios, the overall average NO₃-N load increases, but the magnitudes of the May and November loads are damped. The increase in NO₃-N load in spring appears to occur earlier in the season, while the fall peak occurs later into December. The December peak is higher than the fall peak in the future scenarios; however, these are reversed in the historical period.

Agricultural practices were altered to represent conditions in a future climate by extending the growing season and increasing fertilizer application amounts by 10% and 20% for the mid-century and end-of-century periods respectfully. Model outputs indicated only small increases in NO₃-N load to the lagoon as a result of these changes. For the mid-century period, NO₃-N load was increased by 3%, and for the end of century period, it was increased 1% from the original climate change simulation. An additional scenario was simulated where 50% of crop land was converted to coniferous forest (Supplementary Table S9). Even with this drastic change in land use, NO₃-N load was only reduced by 20% in the historical period and by 12% for 2040 to 2100 when compared to the baseline climate change scenario (Figure 5). There results partially reflect the challenges in managing the ‘legacy’ effects of nitrate as aquifers can continue to pollute surface waters long after land-use practices change at the land surface (Tesoriero et al., 2013).

Model simulations incorporating BMPs revealed a small change (0%–8%) in NO₃-N loading when a 60 m buffer and delayed tillage were applied as conservation tools (Figure 5). Some studies in other watersheds have shown more pronounced benefits from the tested BMPs. For example, studies have identified grassed buffer strips as effective tools for reducing N export in surface runoff (Dunn et al., 2011; Heathwaite et al., 1998). Mankin et al. (2007) noted a strong reduction in total N concentration in field runoff and observed that infiltration accounted for >90% of total N removal. It is expected that grassed buffers have more impact on sediment, phosphorus, and pesticide transport where the main mode of transport is through surface runoff and soil erosion (Al-Wadaey et al., 2012; Sweeney and Newbold, 2014) rather than the groundwater-dominated transport in the Basin Head watershed. NO₃-N is soluble and once infiltrated into groundwater is generally unable to be reduced using BMPs that target surface water and soil erosion.

The SWAT+ model findings for Basin Head differ from results for other PEI studies. For example, delayed tillage of forage crops from fall to spring under a 3-year potato rotation on PEI has been shown to reduce N availability in soil and related leaching into the subsurface (Stuart, 2017), without negative effects on soil health or potato yield (Carter et al., 1998). This practice aims to reduce mineralization of forage crops through the winter months to allow for more N carry-over into the growing season. However, in the modeled scenario, this practice should be paired with a reduction in fertilizer application to impact NO₃-N leaching. Sanderson et al. (1999) found that spring tillage improved potato tuber yields even at lower fertilization rates than late fall tillage. Model results did not indicate that this practice would significantly reduce NO₃-N loading over future decades, likely due to the large input of fertilizer still needed to cultivate potatoes and poor N utilization efficiency of the crop. The only appreciable reduction to NO₃-N loading in SWAT+ came from changing the crop rotation to PSB from PBC. Under the historical period, NO₃-N load is reduced by approximately 33% with this BMP; however, during the mid-century and end of century periods, this reduction is lowered to 26% and 27%, respectively (Figure 5). This indicates that crop rotation change for this watershed may yield reduced benefits as climate change intensifies. A future climate may experience competing phenomena where precipitation and hence hydrologic loads increase, yet high temperatures and low soil moisture decrease nutrient loading. These results are consistent with Liang et al. (2019) which found N surplus in soils to be significantly higher under a PBC rotation than a PSB rotation, hence leading to higher levels of N leaching. The N fixing properties of legumes in the PSB scenario allowed for a significant reduction in fertilizer application (50%), which is likely the cause of the reduction in NO₃-N leaching. It should be noted that this is a simplistic approach that does not consider exact N fixing characteristics of crop rotations at the site, and no in situ soil data was used to confirm the application rates.

A main limitation of the model results are that the findings cannot be generalized to watersheds outside the province of PEI due to the unique soil, crop and aquifer characteristics. Also, the climate change projections (GCMs) used were chosen solely based on cumulative precipitation, representing low, mid-range and high cumulative precipitation, thus, may not represent the full range of potential future conditions. Another limitation that should be noted is that the surface water watershed modeled map not necessarily align with the watershed of the aquifer. This is an assumption that is commonly employed in such studies (Fan, 2019). Based on validation efforts using additional streamflow records, high uncertainty exists in all modelled values currently produced. This model with a relatively short calibration period is used to simulate long-term climate impact, which enhances measurement uncertainty which affects the model confidence. Based on these limitations, such a model should primarily be used for exploratory purposes (Harmel et al., 2014), as in the present study. Streamflow rating curves may not account for high flow conditions properly due to most in situ measurements being collected at low flow (Supplementary Figure S4). Thus, extreme weather events cannot be accurately simulated. This introduces uncertainty in NO₃-N and TSS calibration, which uses streamflow to convert from load to concentration (Supplementary Figures S6, S7). Future work to reduce model uncertainty should include the collection of more flow with a higher range of hydrologic conditions, NO₃-N, and TSS data and watershed properties for model calibration and validation. Furthermore, it should also be emphasized that there is uncertainty in the prediction of NO₃-N transport because of the unknown secondary effects of climate change such as land use change, and changes to agricultural practices (i.e., fertilization rates and timing, growing season) (Howden et al., 2007; Jeppesen et al., 2011). For this reason, future studies are warranted that consider a broader range of future agricultural scenarios and seasonal shifts in hydrologic forcing.