Neural network model predictions for phosphorus management strategies on tile-drained organic soils

Organic soils, or Histosols, are used for intensive vegetable production. These soils require high-fertilizer applications due to their inherent low nutrient status (Zheng et al. 2014). However, excess fertilizer leaches into the environment, becoming a major contributor to eutrophication in water bodies (Krasa et al. 2019). Studies have found that a reduction of phosphorus (P) in freshwater systems, more than nitrogen (N), is critical to control eutrophication (Schindler et al. 2008). Instances of excessive P loads from cultivated organic soils have been observed in New York (Longabucco & Rafferty 1989), Florida (Noe et al. 2001) and Ontario (Winter et al. 2007).

The Holland Marsh of Ontario, Canada is an area of intensively cultivated organic soil that was partially converted from a natural wetland for agricultural use in 1925. The Holland Marsh represents 1% of the Lake Simcoe watershed but contributes up to 5% of the total P (TP) load entering the lake and is therefore considered a source of anthropogenic pollution (Winter et al. 2007). Due to the high water holding capacity and higher water tables encountered on organic soils, subsurface tile drains are installed to support crop production (Hallema et al. 2015). Though there have been studies on soil nutrients in Holland Marsh (Zheng et al. 2015), research on drainage water quality is limited (Grenon et al. 2021).

Drainage water management in the form of controlled drainage (CD) has the potential to mitigate the release of nutrients by controlling the water table within a field, and therefore reduce the amount of water released. The CD functions by placing a control structure on the collector tile line within an agricultural field, close to the outlet, which allows for the control of the water table level within the field. CD has been mostly implemented on mineral soils (Williams et al. 2015a), where multiple studies have concluded that CD is a functional management system that reduces drainage and limits nutrient loads (Sanchez Valero et al. 2007; Skaggs et al. 2012; Carstensen et al. 2016). However, there have been no studies on the use of CD on temperate agricultural organic soils.

Field studies of water quality management practices involve time-intensive and costly processes of instrumenting fields and measurement of water quality. Deterministic and statistical models can reduce the costs of data collection and laboratory analysis of water samples and can be used to fill gaps in data and make hydrologic and water quality predictions under various environmental scenarios (Tiyasha et al. 2021). Neural network (NN) models can be trained to predict water quality using various hydrologic and environmental input parameters (Chang et al. 2016; Liang et al. 2020), however, there are limited studies that use P as an input/output (Tiyasha et al. 2020) and even fewer that look at the management of P (Welikhe et al. 2021). NN has an advantage over process-based models (Qi et al. 2018), which require multiple specific ancillary parametric data points to create a robust simulation. NN is a method of data analysis that uses machine learning (ML) and analytical model building to understand the relationship between multiple inputs with the output variable, without detailed knowledge of the internal process-based functions (Ha & Stenstrom 2003). This type of model has been used extensively for basin rainfall–runoff processes (Dahamsheh & Aksoy 2009), nutrient functions in lakes (Sinshaw et al. 2019), and has been used to predict P loads (Chang et al. 2016) in groundwater, streams and rivers. However, the use of NN in predicting nutrient loading in agricultural drainage runoff is limited. Furthermore, optimization algorithms are effective tools for developing cost-effective watershed designs that meet specified water quality objectives (Seppelt 1999). Environmental simulation modeling with an optimization application can include established objectives, such as pollution load reduction, into specific management plans (Wu et al. 2018; Qi et al. 2020).

The study aims to investigate NN model approaches for predicting P loads from the organic soils of the Holland Marsh. Beneficial management scenarios were assessed using the highest performing NN models, combined with an optimization algorithm, for reducing P loads under CD at a field site in the Holland Marsh.

The study was conducted across 2 years (2015–2016) on a 4.2 ha field (44.064900 Latitude, −79.587878 Longitude) with subsurface tile drains within the Holland Marsh, Ontario. Environmental, agronomic, management and climate data were collected over the course of the study and used to calibrate and validate NN models. The NN models were further used to predict water management scenarios to identify beneficial management practices (BMPs) farmers can use to reduce P pollution. This section includes the following: (1) description of the study area, instrumentation, data collection and water quality analysis; (2) NN modeling methodology, parameters, description of NN models, optimization of the models and performance metrics and (3) the direct and inverse approach and creation of the management scenarios.

The study area has a temperate climate with approximately 5 months (November–April) of the year under snow cover. During the growing season (May–October) average temperatures range from 10 to 22 °C (Government of Canada 2018). The monthly precipitation for 2015, 2016 and a 30-year average is presented in Supplementary Table S1. The annual precipitation was 590 and 552 mm for 2015 and 2016, respectively, compared to the 30-year average of 780 mm, suggesting both years were drier than average. The general soil properties for the study site were analyzed by AgriDirect (Longueuil, Quebec) using a composite soil sample taken post-harvest 2015. The soil properties consisted of an organic matter content of 68%, a bulk density of 0.31 g cm⁻³, a pH of 6.6, and a soil Mehlich P of 655 mg kg⁻¹. Carrots (Daucus carota Bergen) were cultivated for the duration of the study. The mineral fertilizer application in May 2015 consisted of 20 kg-P ha⁻¹ and 210 kg potassium (K) ha⁻¹, while in May 2016 it consisted of 10 kg-N ha⁻¹, 20 kg-P ha⁻¹ and 205 kg-K ha⁻¹. Seeding occurred soon after the fertilization, and the crop harvest started in September for both years.

A CD structure (Supplementary Figure S1) was installed at the outlet of the collector tile line. This structure allowed for manual control of the water table height in the field with the placement of removable gates that stack on top of each. The water level strategy adopted was to have the water table between 30 and 40 cm from the soil surface during the non-growing season (November–April) to reduce tile flow during the spring, and between 70 and 80 cm from the soil surface during the growing season for agronomic management. Similar strategies were implemented in other studies where drainage reduction during the spring-thaw period is critical for nutrient load reduction (Gunn et al. 2015). The CD gates were removed in May to allow water to drain from the field for fertilizer application and seeding. In 2015, the height of the gates varied from 70 cm to the drain depth of 1.48 m from the soil surface between mid-June and early July to reduce surface water ponding due to excess rain. While in 2016, they were lowered mid-April to dry the field for planting and then raised in July to conserve water within the field due to lack of precipitation.

Discharge measurements were taken using a compound weir in the CD structure, which was placed as the uppermost gate. The compound weir consisted of an 11° V-notch to estimate flows less than 1 L s⁻¹ and a rectangular weir to calculate flows above 1 L s⁻¹. The discharge was continuously recorded at 15-min intervals using a pressure transducer within the control structure.

Water samples were collected during periods of drain discharge, using an ISCO 6712 portable auto-sampler at 4-h intervals (Teledyne ISCO, Lincoln, Nebraska). Weekly discrete grab samples were also taken. No samples were taken during periods of winter freeze or periods of no drain flow due to dry weather conditions. The water samples were stored in a 4 °C refrigerator at the Guelph University Muck Crops Research Station (MCRS) before being transported in coolers with ice to the McGill University laboratories for analysis.

Water analyses consisted of TP using the potassium persulfate digestion method (Dayton et al. 2017). For each analysis, both blank distilled water samples and standards of P were included for the calibration of the results. Calculation of P loads was done using TP concentrations (mg L⁻¹) and continuous discharge (L s⁻¹). The load was calculated using linear interpolation, which is one of the most common methods of load calculation with the least amount of uncertainty (Williams et al. 2015b).

The non-linear or stationary feature of water in addition to its interdependent relationship with the environment and agronomic interactions creates a complex problem (Tiyasha et al. 2020). Qi et al. (2018) reviewed biophysical phosphorus models and found that the performance was dependent on a complex assessment of environmental parameters. Furthermore, these models have been calibrated for mineral soils and there is a lack of information on their performance with organic soil (Grenon et al. 2021). However, NN models allow for the incorporation of relevant and available environmental data without having to account for the entire physical environment and its processes. NN models have been extensively studied and concluded as a reliable tool to deal with non-linear water quality data, to increase decision-making capacities and to study various water quality management problems (Tiyasha et al. 2020).

An ad hoc model-free approach (Maier et al. 2010) based on domain knowledge was used to select input parameters. The daily input parameters included maximum temperature, minimum temperature, precipitation, evapotranspiration, soil moisture, the water level within the field, water height within the CD structure, CD discharge and fertilizer application. A comparison analysis was conducted to assess the interaction of the parameters with the model output, TP loads.

Seven NN models were chosen for this study, using two modeling approaches. The seven models consisted of feed-forward NN (FNN), deep feed-forward NN (DFNN), long-short-term memory (LSTM), bi-directional LSTM (Bi-LSTM), gated recurrent unit (GRU), general regression NN (GRNN) and radial basis function NN (RBFNN). All models were run using two different modeling approaches: a time-series and a randomization analysis. The time-series analysis was done using 2015 data to train the model and 2016 to test the accuracy of the model. This allowed for the inclusion of temporal trends in the analysis. The randomization analysis was done by randomly assigning 80% of the 2 years of values to train the model, and the other 20% of the values to test model accuracy. In total, the 14 (7 time series and 7 randomizations) models were run on a daily time scale.

A feed-forward neural network (FNN) model includes one hidden layer, while multiple hidden layers are called a DFNN. The first and the last layers are for inputs and outputs, respectively. The operating principle of a typical FNN is based on the behavior of neurons and synaptic connections, where the input data are processed through three different operations: (i) multiplication by the synaptic connection weights; (ii) summation by adding the threshold function and (iii) transformation of results by applying an activation function. In the training phase, FNN automatically adjusts their synaptic weights and bias through back-propagation (Schmidt et al. 2018). The Bayesian optimization procedure was to select the optimal training function. Supplementary Figure S2 shows the schematic representation of typical FNN and DFNN structures while the equation in Supplementary Table S2 characterizes the relationship of one neuron.

Recurrent neural networks (RNNs) are based on the concept of time layering, which gives them the ability to take temporal series data as inputs and outputs. In more general terms, RNN is particularly suitable for sequence-to-sequence learning purposes (Sainath et al. 2015). The standard structure of an RNN model consists of an input layer, multiple hidden layers, and an output layer. In contrast to other RNN, memory cells, which are the main features of LSTM and Bi-LSTM networks, are incorporated into the hidden layers. The equations governing memory cell operations can be found in Supplementary Table S2. Memory cells are composed of an input (i_t) gate; to determine which information will be added to the cell, an output (o_t) gate; to identify which information can be used as output in the cell and a forget (f_t) gate; to decide which information can be removed from the cell. At each time step, the input and output are introduced into the three gates and candidate values are computed (Supplementary Figure S3).

In a similar manner to the training phase in FNN, the weights and bias terms are computed using the back-propagation algorithm with gradient descent so that an objective function is minimized across temporal target values. The idea behind the Bi-LSTM is to aggregate input signals in the past and future of a specific time step in LSTM models. The GRU is an RNN similar to LSTM. However, in contrast to LSTM, GRU has only two gates: the update and reset gates which decides what information is relevant and should be passed to the output.

These ML models are neural networks with radial basis functions as activation functions. Global approximation, compact scheme, and ability to fit continuous and noisy data are features that distinguish RBFNNs and GRNNs from other NNs (Yu et al. 2011). The most used RBFNN activation function is the Gaussian function, which possesses only a few hyper-parameters. RBFNN's principle of operation does not differ from FNN. However, the training process is tuned by incorporating some clustering techniques such as decision trees (Kubat 1998), or k-means (Sing et al. 2003). The typical structure of RBFNN includes a two-layer network: the first layer comprises RBF neuron cells that compute weighted inputs and process them to the second layer with a linear transfer function. The only difference between RBFNN and GRNN is the use of biases in both layers of RBFNN and no biases are used in the case of GRNN (Ozyildirim & Avci 2013).

Optimization is conducted within the model to acquire the best possible results between the training of the model and the testing of the model accuracy. The range of parameters and hyper-parameters used to build the NN models are presented in Supplementary Table S3. Hyper-parameters are vital to control parameters for ML where their appropriate choice leads to high-performance training and thus highly effective predictions. Many optimization techniques have been developed to tune ML hyper-parameters. A Bayesian optimization approach was the most suitable, as a probabilistic model with prior distribution is built through past evaluation. This meta-model surrogates the objective function for a much easier optimization (Snoek et al. 2015). The model optimization occurs through iterations. Each iteration allows for different configurations between the learning rate, and the number of epochs within the FNN, DFNN and RNN models. The GRNN and RBFNN model iterations allowed for the spread of the values to be varied and optimized. All models were run at least three times using 100 iterations to gain the best optimization possible.

The performance of the models was evaluated using the following statistical measures: Nash–Sutcliffe efficiency (NSE), correlation coefficient (R) and root mean squared error (RMSE). Furthermore, the comparison of the measured and the predicted TP loads were evaluated with the percentage of uncertainty to identify the accuracy of the models (Williams et al. 2015b).

Through the optimization of the various models, the best model was found for each approach: one for time-series and one for randomization. Although all models were optimized, predicting TP management strategies was performed using the best model for each approach. Two TP management strategies were used: a direct approach and an inverse approach. All management scenarios are presented in Table 1. The management scenarios were run by separating the spring (January–April) and summer (May–October) periods. This was done because the water management strategies differed between these periods due to agronomic activities. The objective was to find the BMP for water management, aimed at reducing TP.

The direct approach (Supplementary Figure S4) was done using the calibrated models and running new input parameters through the models to predict the TP load. The new input parameters used the nine, 2015 input datasets and varied the values within the file depending on the management scenario (Table 1). The spring and summer management scenarios were done by re-calculating what the discharge would be at different set gate levels within the CD structure. As the discharge is directly related to the water table height, this allows for a comprehensive understanding of how the management practices affect the total nutrient load. The six spring scenarios were run by varying the set gate-level between 20 and 45 cm from the soil surface, while five summer scenarios were run, maintaining a constant water height while accounting for the 9 days of free drainage needed to reduce ponding due to rainfall. Other management scenarios included four fertilizer rates, from no fertilizer to five times the actual rate, to simulate the effect of fertilizer application on TP loads. This consisted of retaining all 2015 input parameters except the fertilizer rate which was manually changed to reflect the new application rate. Finally, using the optimal management strategy found from the spring and summer scenarios, these new set gate heights were applied to a year with 900 mm of rain (wet year: 15% greater rainfall than the 30-year average), and the 2015 data which was a dry year (Supplementary Table S1). Climate data from 2010 were used as the input parameters for the wet year.

The inverse approach with particle swarm optimization (PSO) was the second management strategy used. The PSO's mechanism was inspired by the swarming behavior of biological populations (Eberhart & Shi 2001). The computational procedure of the PSO algorithm is described in the following steps: (i) initialization of the set of PSO particles with random positions and velocities; (ii) evaluation of the objective function for each particle; (iii) the value of TP load is calculated at each new particle position until a better position is achieved by the PSO particle and the value of P_best (best position) is replaced by the current value and (iv) testing the model that if the new value of the objective function (called fitness) is better than the old, the value is updated and stored in a list. In brief, this inverse method allows for the PSO algorithm to be coupled to the NN model and predict an input parameter while maintaining a specifically wanted output.

The proposed discharge by the PSO is applied as an input variable with the other eight 2015 input variables to create a new input dataset. This input dataset was then entered into the optimized NN model to produce the corresponding TP load. The resulting TP load is compared to the calculated TP load. If they are equal, then this computed TP load is stored as a potential solution to the inverse analysis problem (Supplementary Figure S5). In contrast, if the values of the desired and computed TP load are different, the PSO iterations are run until convergence.

The PSO approach was used with various management scenarios (Table 1) which consisted of reducing the TP loads found in 2015 by a percentage varying between 10 and 60%. These new values were then applied to predict the discharge and set the gate-level within the CD structure needed to achieve this reduction.

The TP load (Supplementary Figure S6) totaled 0.45 kg ha⁻¹ in 2015 and 0.50 kg ha⁻¹ in 2016. Tile drain discharge occurred during spring-thaw each year and high rainfall events. Studies have shown that some of the most significant discharge events occur during spring-thaw (Lam et al. 2016). This period of high discharge was corroborated by the 2016 results (spring load: 0.33 kg ha⁻¹; summer load: 0.17 kg ha⁻¹), where 67% of the tile discharge occurred from February to April resulting in 66% of the TP load, due to the freeze-thaw, as well as the increased precipitation in March 2016 (80 mm), compared to the 30-year average (46 mm). Furthermore, the 2015 TP load (spring load: 0.18 kg ha⁻¹; summer load: 0.27 kg ha⁻¹) was influenced by fertilizer inputs followed by high rainfall in June 2015 (160 mm), which was 112% higher than the 30-year average (75 mm), resulting in 60% of the TP load occurring during the summer. The annual TP load in the drainage water accounted for approximately 3% of the 20 kg-P ha⁻¹ applied fertilizer for both 2015 and 2016. Therefore, only a small portion of the applied P is released through tile drain discharge (Sanchez Valero et al. 2007). Large amounts of P were sequestered within the soil or taken up by crops, as found in other studies (Easton & Petrovic 2004).

There was a high correlation (0.97) between the discharge and TP load (Table 2). As the TP load was calculated using both the continuous tile discharge and the TP concentration (Williams et al. 2015b), this high correlation was expected. A correlation was further found between the climate parameters (temperature, evapotranspiration) and the soil hydraulic parameters (soil moisture, field water table). Other studies have found that temperature and precipitation can impact P movement (Liu et al. 2021). The relationship between temperature, evapotranspiration and soil properties is reflected in the model. Soil hydraulic parameters and water levels were not significantly correlated with precipitation. Although climate parameters are important to understanding nutrient movement, trends are not often discernible (Coffey et al. 2019). Sinha & Michalak (2016) concluded that spatial variability of nutrient loading is driven by nutrient inputs, but intense precipitation events, greater than the average rainfall, are a driver to increased loading. The accuracy of NN models can be improved when there is an understanding of the functional relationship between the input parameters (Tiyasha et al. 2020). Therefore, the input parameters were chosen with an understanding of their inherent relationship to the output TP load.

The time-series and randomization analysis for FNN and DFNN (Table 3) resulted in higher accuracy under the latter analysis. Others have found that FNN models have better performance compared to other NN models (Ahmed 2017). The simple structure of one hidden layer optimizes the linear parameter interactions and creates the optimal solution (Tiyasha et al. 2020). Studies have found that DFNN is an improvement over FNN models (Wilamowski & Yu 2010). As shown in Table 3 for randomization, the performance of the model increased between FNN (R=0.92) and DFNN (R=0.98). However, the uncertainty also increased under DFNN (8.54%). Training of DFNN on a small dataset caused increased uncertainty (overfitting problem) within the model prediction. Studies have found that the overfitting in training DFNN models can cause an accumulation of errors within the neurons, and decrease the model accuracy (Mhaskar & Poggio 2016).

The time-series models (Table 3) resulted in low R values (FNN: 0.58; DFNN: 0.48) and NSE (FNN: 0.5; DFNN: 0.38) and high percentages of uncertainty (FNN: 63%; DFNN: 59%). The inaccuracy of the time series can be accounted for by the type of NN model. The FNN models need more data for training (Chen et al. 2017). The time-series analysis, however, had a higher percentage of data points used for model testing, compared to the randomization analysis.

The RNN models did not perform well using a randomization analysis. The RNN models can process temporal data identifying interconnections and allowing feedback more than other NN models (Chen et al. 2020). Therefore, the temporal aspect of the time-series analysis is necessary for RNN models. The time-series analysis (Table 3) showed that the LSTM model outperformed both the Bi-LSTM and the GRU models. In addition to having better model performance in comparison to Bi-LSTM and GRU, the LSTM model had an overall lower percentage of uncertainty at 1.64% (Table 3). Other studies have found that LSTM has a good performance when analyzing nutrient parameters in water and air quality (Hamrani et al. 2020; Huan et al. 2021). Studies have found that the added ‘memory cell state’ within the LSTM model allows for a greater understanding of the input data and higher accuracy (Wang et al. 2020).

The results obtained with the RBFNN and GRNN models showed both underperformed with the time-series (GRNN, R=−0.45; BRFNN, R=0) and randomization (no results) analyses. Essentially RBFNN and GRNN models are very sensitive to the high variation of input and output data, as well as to the data distribution (Heydari et al. 2019). These models are, therefore, more suitable for normally distributed data. Furthermore, the low accuracy of these models can be linked to poor flexibility in terms of tunable hyper-parameters, as a single parameter can be optimized (Benoudjit & Verleysen 2003), unlike the FNN or RNN models that have multiple optimization options. Furthermore, the dataset of 2 years can be considered as the minimum data required for model predictions and could, therefore, have impacted the performance of these models as past studies have found RBFNN an accurate model to use in water quality predictions (Sharma et al. 2003).

The FNN was found to be the most accurate model under a randomization analysis, while under time-series, the LSTM outperformed the other models. The results from the NN models identified that the overall performance of the optimized FNN and LSTM models was high and allowed for the prediction of management scenarios. Two management strategies were undertaken: a direct approach and an inverse approach. The management scenarios (Table 1) allow for an understanding of BMPs and how varying these practices can affect the TP loads in the environment. The management scenarios were run to predict the most effective water table height during the spring and summer periods to reduce TP loads.

The direct approach for predicting TP loads found that the LSTM predictions had an uncertainty of 13.2% and the FNN 6.2% compared to the actual TP loads of 2015 (Figure 1). Furthermore, the water table management scenarios (Table 4) were overall more accurate for spring-thaw compared to the summer period where all predictions were significantly lower than actual TP loads. Chen et al. (2020) found that at minimum 1 year of data is needed for long-term forecasting. Therefore, this study had available the minimum data required. The predicted peaks of TP loads (Figure 1) from the FNN model were similar to the measured, though not as high a daily value.

Table 4 suggests that setting the water table at 35 cm or lower from the soil surface can increase TP loads compared to a water table of 30 cm. The management of water levels during the spring is crucial, as most of the nutrient discharge occurs during this period (Williams et al. 2015a). Therefore, maintaining a higher water level can reduce the annual nutrient loads. The spring FNN results were more consistent than the summer, where there was less perceived change in TP loads between 80 and 50 cm. The summer LSTM water table management scenarios showed that 70 cm from the soil surface can reduce the TP loads by 32%, compared to 80 cm or deeper. McDonald & Chaput (2010) noted that carrots require the water table to be between 75 and 90 cm for optimum growth. Therefore, given the difference in hydraulic gradient between the water level in the CD structure and the upstream tile lines, our predictions for a water table height of 80 cm from the soil surface in the CD structure are well in line with the recommendations of McDonald & Chaput (2010).

The results (Table 4) found no variation in TP loads between actual fertilizer rate and no application of fertilizer. As all other climate and hydrological parameters were the same, the models are unable to predict fertilizer effects when the rate is decreased. The results identify increased TP loads with the FNN model at double the fertilizer rate with a range of 0.43–0.57 kg ha⁻¹ and at five times the fertilizer rate, with a range of 0.57–0.89 kg ha⁻¹. Studies have shown that nutrient inputs are the driving factor to nutrient load increase (Sinha & Michalak 2016), which is seen in increased loading with increased fertilizer application.

Table 4 shows that the BMP under the direct approach is a spring water table set at 30 cm from the soil surface, and 80 cm during the summer period. The FNN model predicted a TP load of approximately 0.21–0.22 kg ha⁻¹ during a dry year with a slight increase in the rate (0.23–0.41 kg ha⁻¹) during a wet year. The LSTM predictions for the dry (0.24–0.30 kg ha⁻¹) and wet years (0.35–0.46 kg ha⁻¹) were higher compared to the FNN. There is a 15% increase in rainfall between the dry and wet years, signifying an increase in nutrient loads. However, the predictions identify similar loads between the wet and the actual 2015 practice. This identifies that free drainage and precipitation, even for a short period (as seen in 2015), increased nutrient load. Studies have found that RNN models can better handle historical time-series input data to predict water quality (Li et al. 2019). Therefore, according to past studies, the LSTM predictions may be more accurate.

The inverse approach with the PSO algorithm quantified the TP load reduction rates and estimated the optimal corresponding discharge. The optimum values of water discharge corresponding to the TP reduction rates (10–60%) obtained from FNN-PSO and LSTM-PSO models during the summer season are presented in Figure 2. The predicted discharges were not able to identify the peaks caused by excessive rain in June/July 2015. The mean predicted discharge under TP load reduction (%) and the subsequently calculated gate levels are given in Table 5. In the spring season, the reduction in average discharge is higher with the FNN models, however, the percent reduction of TP was similar for each water table height. A reduction in TP load of 20% showed similar results under the direct approach (Table 4) reflecting that the optimal water table level from the soil surface is approximately 30 cm. Water quality management methodologies based on the PSO approach have been successfully applied for optimizing BMPs (Liu et al. 2020; Rohmat et al. 2021). However, only a few studies used the combination of ML models with optimization approaches (He et al. 2020; Taghizadeh et al. 2021). This inverse analysis resulted in the calculation of the necessary water table height for the desired TP load and can therefore facilitate the adoption of BMPs at a field scale.

The summer period results (Table 5), with the FNN-PSO analysis, showed that decreasing the TP load by 10% is obtained by a reduction of 35% in the average drainage discharge. This reduction in discharge, however, was not found with LSTM-PSO. Both FNN-PSO and LSTM-PSO found that the TP load reduction of 20–60% identified a gradual decline in discharge. The PSO results predicted that reducing the TP load by 10% would set the water table level within the CD structure to approximately 85 cm from the soil surface, while with a reduction of 60% TP loads, the water table level would be about 71 cm. As carrot crops need a water table level between 75 and 90 cm from the soil surface (McDonald & Chaput 2010), the inverse analysis predicted that a reduction of 30–40% TP loads would still be obtained with a water table depth of approximately 78 cm and would be able to meet the crop aeration requirement by adequate drainage.

This study investigated the optimization of seven NN models (FNN, DFNN, LSTM, Bi-LSTM, GRU, RBFNN, GRNN) under two different modeling approaches (randomization and time-series). Our results show that the FNN and LSTM models are very reliable with only a few years of field data. Continued analysis into parameter distribution and longer-term datasets can obviously strengthen model performance. The results identified LSTM as the model with the best performance under the time-series approach with temporal time step being important to the performance of RNN models. In addition, the FNN model outperformed all other models under the randomization approach. Two management strategies, the direct and the inverse approach, were used to predict the best water management strategy needed to reduce TP loads. Both the direct and inverse approach predicted similar recommendations that can be adopted by farmers. The direct approach found that the optimal strategy was to set the water table at 30 cm from the soil surface during the spring-thaw period and at 80 cm during the summer period, to reduce the TP load to the environment, while also not negatively affecting crop production. The direct approach was also used to assess fertilizer management scenarios and found that any increase in the fertilizer rate from actual (20 kg-P ha⁻¹), will increase the TP load. The second management strategy, the inverse approach found that a reduction of 20% in the spring TP load can be achieved by raising the gate-level to 30 cm from the soil surface. Furthermore, in the summer period, the analysis identified that a reduction of up to 40% of the TP loads was possible while maintaining a water table height of approximately 78 cm for crop production.

The broader applicability of our research lies in the ability to apply the LSTM or FNN models to multiple locations within an organic soil agricultural area, thereby allowing for the development of appropriate water management and fertilizer application strategies to achieve P reduction in lakes and rivers. NN models can optimize fertilizer rates and water table levels for different crops and fields, based on inputs of climate and soil moisture data for P management strategies.

We thank Hicham Benslim, Hélène Lalande, David Meek and Khosro Mousavi for their technical assistance and are grateful to Aidan De Sena, Paddy Enright, Naresh Gaj, Samuel Ihuoma, Marjorie Macdonald, Kerri Edwards, Kenton Ollivierre, Naeem Abbasi and Rebecca Seltzer for their field/lab assistance. We also thank the University of Guelph-Muck Crops Research Station staff and the landowners who provided access to the fields during this study.

Funding for this project was provided by the Natural Sciences and Engineering Research Council of Canada (NSERC) under the Strategic Projects Grant fund # 447528-13.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.