Simulated effects of gypsum amendment on phosphorus losses from agricultural soils

Elina Jaakkola1*, Sirkka Tattari1, Petri Ekholm1, Liisa Pietola2,Maximilian Posch3 and Ilona Bärlund4 1 Finnish Environment Institute, P.O. Box 140, FI-00251 Helsinki, Finland 2 Yara International ASA / Yara Suomi Oy, Mechelininkatu 1 a, FI-00181 Helsinki, Finland (present address: Central Union of Agricultural Producers and Forest Owners, P.O. Box 510, FI00101 Helsinki, Finland) 3 Coordination Centre for Effects, RIVM, P.O. Box 1, NL-3720 BA Bilthoven, The Netherlands 4 Helmholtz Centre for Environmental Research-UFZ, Brückstrasse 3a, D-39114 Magdeburg, Germany * e-mail: elina.jaakkola@environment.fi


Introduction
Phosphorus (P) is a necessary nutrient in crop and animal production.Because of its losses from fields to surface waters, agriculture is a major contributor to eutrophication (Sharpley andRekolainen 1997, Carpenter et al. 1998).Despite several measures to control P losses (e.g., reduced tillage methods, lower fertilisation rates, vegetative buffer zones, and constructed wetlands), reduction of P losses has been slow, and, therefore, novel methods are sought.
To the best of our knowledge, the effect of gypsum has not been incorporated into nutrient transport simulation models.The field-scale model ICECREAM (Rekolainen andPosch 1993, Tattari et al. 2001;Yli-Halla et al. 2005) has been used to predict the effect of other agri-environmental measures, such as fertilisation, tillage practices, and buffer zones, on P losses in Finland (Rekolainen et al. 1999, Tattari et al. 2001, Bärlund et al. 2009) and in Sweden (Arheimer et al. 2004).The objective of this study was to add the effect of gypsum to the ICECREAM model.The amended model was calibrated and verified against data from a small agricultural catchment in southern Finland, where gypsum was spread and its effects were monitored by means of automatic water-quality sensors and grab samples.In the previous version of ICECREAM (ICE0608), water was routed via surface runoff and through the soil profile as matrix flow.However, Finnish clay soils with subsurface drainage show short response times between rain events and drain flow as well as high nutrient concentrations in the drain flow soon after fertilisation, which suggests that macropores contribute to drain flow (Paasonen-Kivekäs et al. 2008).Models such as CROPWATIN and MACRO, when applied to Finnish clay soil, estimated the proportion of macropore flow in total subsurface drain flow to be as great as 87-99% for fields where subsurface drain flow's share in the total drain flow was measured to be 55-64% (Karvonen 2004, Hintikka et al. 2008).
In Sweden too, macropore flow has been found to be a major flow path and a description of macropore flow has already been included in the locally tuned version of the ICECREAM model (Larsson et al. 2007).However, the Swedish macropore model introduced a new groundwater reservoir separated into two compartments, which requires data (e.g., the amount of readily available dispersible particles and soil detachability coefficient) that have not been measured in Finland.Therefore, a simpler approach was introduced for our purposes through modification of the variables and parameters for factors affecting erosion and P transport.
In this paper, the ICECREAM model, extended to include macropore flow, was first calibrated for a situation before gypsum amendment.After this, the effect of gypsum was added to the model on the basis of laboratory experiments performed in the TraP project (e.g., Uusitalo et al. 2012) that suggested decreasing P losses and in view of field experiments (Cox et al. 2005, Hämäläinen et al. 2010) suggesting changed hydraulic properties of the soil.Other laboratory experiments have also shown that gypsum affects < 2 µm sized aggregates, probably by decreasing so-called diffuse erosion (see Aura et al. 2006) that may take place when the soil is wet (e.g., when soil frost has melted and released water to the soil) and small clay particles can be transported by Brownian motion from aggregates to soil water.However, a change in the stability of microscopic aggregates is difficult both to measure and to add to any model.Thus, other technical approaches, such as lowering the concentration of soluble P in soil water and changing the hydraulic properties of the soil, were selected to simulate the gypsum amendment.The modified model was verified for the catchment experiment, wherein gypsum was amended and P losses monitored.

Brief description of the ICECREAM model
The ICECREAM model is based on the CREAMS model (Knisel 1980), modified for Finnish conditions (Rekolainen and Posch 1993, Tattari et al. 2001, Yli-Halla et al. 2005, Bärlund et al. 2009).The hydrology is calculated as daily runoff via a modification of the curve number method (Soil Conservation Service 1972).The curve number depends on the soil cover and on the hydrological properties of soil types.First the amount of surface runoff is calculated, after which the rest of the water is assumed to infiltrate.Evapotranspiration is computed via the Penman-Monteith model (Monteith and Unsworth 1995).The erosion sub-model is based on the modified USLE (Foster et al. 1977, Posch andRekolainen 1993).The model is developed for mineral soils and should not be applied to soils with high organic matter content.
The P sub-model is based primarily on the GLEAMS (Knisel 1993) and EPIC (Jones et al. 1984) models.Soil P is divided into different pools, as in GLEAMS, which are labile P, active inorganic P, stable inorganic P, fresh organic P, and more stable organic P. Initial values must be assigned to all pools for each soil layer.The total P (TP) lost from the soil in ICECREAM consists of PP and DP (Tattari et al. 2001), which is not further separated into reactive and non-reactive P. The new features included in the model are described below.

Incorporating macropore flow into ICECREAM
The ICE0608 version (Bärlund et al. 2009) of the ICECREAM model only simulated percolation via matrix flow.The concentration of percolating DP was calculated from the deepest layer of the soil and is a product of the entire soil profile.However, models not accounting for macropore flow may overestimate P losses via surface runoff and erosion in areas where considerable amounts of water are routed via subsurface drains rather than via surface runoff (e.g., Larsson et al. 2007, Hintikka et al. 2008, Warsta et al. 2009).In the ICE0609 version of ICECREAM, developed here, macropores are considered as a simple transport pathway, in which all P processes are lacking.The water flow is first divided into surface and subsurface flows.The subsurface flow consists of matrix flow and macropore flow, with the latter occurring only when the soil is moist enough (see below).Only the two topmost soil layers (the top 10 cm) were assumed to contribute P to macropore flow.
For simulation of the macropore flow, three new parameters were introduced to the model: The relative area of macropores (0-1; e.g., 0 for no macropores and 0.05 if 5% of the soil layer consists of macropore inlets) controls the amount of water routed via macropores (instead of matrix flow).Typically, the area of macropores (diameter > 0.03 mm) decreases with depth, being 2.8-11% in the 0-20 cm layer but only 0.3-2% in the 30-50 cm layer (Aura 1995, Pietola and Yli-Halla 2007, Hintikka et al. 2008, Alakukku et al. 2010).However, in the ICE0609 version, this parameter is the same from top to bottom.Second, a dimensionless threshold value in the two uppermost soil layers (10 cm) indicates the soil moisture level above which macropore flow starts.Finally, a dimensionless ω factor determines the percentage of PP transported from the P pools to macropores.Its value does not depend on other factors, such as the amount of water.
The DP (kg m -2 ) in surface runoff is computed for the uppermost soil layer as C Pil,av = concentration of labile P in the soil water of the uppermost soil layer (g kg -1 ) q surf = daily runoff (mm day -1 ) β = extraction coefficient (-) C Pil = concentration of labile P in the specific soil layer (g kg -1 ) I = daily infiltration α = initial absorption of water, which depends on soil water content (-) ε = porosity of the surface layer (fraction) ρ s = specific density of the surface layer (kg dm -3 ) K d = partitioning coefficient (dm 3 kg -1 ) Percolation of P (kg m -2 ) out of the matrix is computed from the concentration of P in the deepest layer in the soil profile (the same equation applies to the other layers): (3) where (4) and q perc = daily runoff through matrix (mm) P il = labile P in a specific layer (kg m -2 ) m s = soil mass (kg m -2 ) sw = total amount of water in the soil layer (mm) The DP in macropore flow -that is, the bypass (kg m -2 ) -is computed by means of the concentration of P in soil water in the uppermost two layers (10 cm), where the concentration is calculated via equations 2 and 4, with the amount of water (q bypass ) comprising both the soil water and the water via macropores: (5) where q bypass = daily runoff through macropores (mm) The PP in surface runoff (kg m -2 ) is calculated as the sum of the labile P pool (PP il ), active mineral P (PP ia ), stable mineral pool (PP is ), organic humus pool (PP oh ), and manure pool (PP om ): (6) D = soil loss (kg ha -1 ) er = enrichment ratio (-) The PP in macropore flow (again, bypass, as kg m -2 ) is computed from the P pools of the two topmost layers: (7) ω = fraction of PP moving from P pools to macropores

The catchment experiment in Nummenpää
The ICECREAM model was applied to the Nummenpää catchment (total area of 2.45 km 2 , with 41% being fields and the remaining area mainly forest), a sub-basin of the river Lepsämänjoki, in southern Finland (see Fig. 1).The catchment area was delineated from a digital elevation model (DEM, National Land Survey of Finland, Licence 7/MML/11) on the basis of LiDAR data (minimum vertical resolution of 0.3 m, grid cell size of 2×2 m).In 2008, there were 73 inhabitants of the study area (YKR 2011), with none of the houses connected to central sewage systems (Autio K., personal communication).The upper, forested parts of the catchment consist of rocky areas (dystric leptosol) and eskers (haplic podzol), with the soil of the lower, agricultural part of the catchment being fine clay (vertic cambisol, Finnish Soil Database 2010).Soil sampling and farmers' own data indicate that the topsoil texture of the fields is mainly silty clay loam.The soil exhibited cracks when drying after spring snowmelt.Calculated via the DEM, the mean slope for the entire field area is 1.6% (minimum mean slope of a field: 0.6%, maximum: 4.8%).
The field area consisted of 32 individual fields, varying in size (0.1-9.9 ha, mean 3.2 ha) and shape, which were further divided into smaller parcels with regard to agricultural practices.In total, 65% of the fields were autumnploughed, the rest having shallow cultivation (32%) or direct sowing (3%) as the tillage method.According to the field register of the Agency for Rural Affairs, the dominant crops in 2009 were spring cereals (~81%) and cabbage (~15%; see Table 1).Soil P status (soil test P, or STP) was analysed through extraction of P with a solution of 0.5 M ammonium acetate and 0.5 M acetic acid at pH 4.65 (Vuorinen a Mäkitie 1955).The average STP, according to soil sampling and farmers' data, was 16 mg l -1 for spring cereal fields (range: 10.8-75.4mg l -1 ) and 50 mg l -1 for cabbage fields (range: 26.0-86.0mg l -1 ).Both values are above the Finnish average for agricultural land (i.e. 12 mg l -1 ), but cabbage cultivation typically receives high amounts of P in fertilisers.The Field Drainage Association of Southern Finland did not have information on the location, directions, and condition of subsurface drainage for every field, but, with the help of the local farmers, the subsurface drainage coverage was approximated to 72%.Most of the drainpipes were installed in the 1950s and 1960s, but according to the farmers, some renovation has taken place in recent decades.

Gypsum treatment and source apportionment
Gypsum was spread on fields (4100 kg ha -1 ) in autumn 2008 after the harvest but before tillage; in total, 91% of the Nummenpää field area was treated.Runoff quantity was monitored at a V-notch weir and runoff quality via grab samples (e.g., turbidity, TP, DP, and DRP), and via an automatic water-quality sensor (YSI 600 OMS for items such as turbidity) before, during, and after the gypsum amendment.For a more detailed description of the monitoring sites and sampling strategy, see the work of Ekholm et al. (2012).
Since ICECREAM simulations were calibrated and verified against the PP and DP loads from the entire Nummenpää catchment, the contribution of fields to the total load had to be estimated.The figures for load from the nonagricultural land were based on grab sampling results upstream of the field area.The load from scattered settlements was based on per capita wastewater P emissions (Onsite Wastewater System Decree 209/2003) and with assumption of 60% retention for this P (Rontu and Santala 1995).With these assumptions, fields accounted for 79% of total P load before the gypsum treatment (Ekholm et al. 2012).Hydrological variation probably affects the P load caused by wastewaters, as revealed by the fact that the main channel is occasionally dry.Because there was no information on the temporal variation of the wastewater loading and its amount was assumed to be negligible in comparison to agricultural loading, its effect was excluded here.
The turbidity (in NTU) measured by automatic sensors was correlated with PP concentrations analysed from the grab samples (as TP less DP) both before and after the gypsum amendment (n = 119), and hourly PP concentrations were calculated as PP = 1.36•NTU + 7.9 (r 2 = 0.92) (8) The DP concentration analysed from the grab samples depended on the flow (Q) before the gypsum amendment (n = 12), and hourly DP concentrations were calculated as After the gypsum amendment, DP showed no correlation with flow.Therefore, a constant concentration was used, corresponding to the average DP concentration (41.8 µg l −1 ) of the grab samples taken after the gypsum amendment (n = 26).The PP and DP concentrations were then converted to loads, and the load corresponding to nonagricultural land was subtracted from the total load by means of a load estimate obtained from a sampling site in the forested upper reaches (see Ekholm et al. 2012) deemed representative of the non-agricultural land.Hereafter, this load, estimated for agricultural land of the Nummenpää catchment, will be referred to as "monitored load".

Other gypsum experiments
In the TraP project, the effect of gypsum was also examined in the laboratory.DRP in percolation water was studied from soils of the Kisko farm fields in southern Finland, which have a barley cropping history of nearly a decade of minimum tillage.Soils were sampled from seven parcels, of differing P status, from 0-10 cm depth after harvest.One sub-sample of each soil (0.5 l) was mixed with 2 g of gypsum (moisture 18% on a mass basis, with the gypsum proportion being comparable to the field applications in the Nummenpää catchment -i.e.4100 kg ha - 1 ), and the other sub-sample was kept as a control.Thereafter, the samples (three replicates for each treatment) were kept close to saturation in plastic pots of 10 cm height, perforated at the bottom.After two or five days' incubation, the soils were watered with 150 ml and the DRP concentration of the percolated water was measured after filtration (Whatman nuclepore polycarbonate, pore size 0.4 µm).
In addition to the soils of the Kisko farm fields, other gypsum experiments were utilised to add the effect of gypsum to the ICECREAM model.Uusitalo et al. (2012) studied erosion and P mobility after gypsum application (3000 or 6000 kg ha -1 ) on two clayey fields (pH of 6.5, STP 10-15 mg l -1 ) that were autumn-tilled to a depth of 20 cm or 10 cm.The study was done with soil cores that were subjected to laboratory rainfall simulations.Hämäläinen et al. (2010) studied the aggregate stability and infiltration capacity under field conditions on fine-textured clay soils (pH 6, STP 9.6 mg l -1 ) that were treated with 0, 2000, 4000, and 6000 kg of gypsum per hectare.The fields were under either direct drill or minimum tillage.Cox et al. (2005) determined the effect of a large-quantity (15000 kg ha -1 ) gypsum application on the timing, quantity, forms, and pathways of flow of P from a grazed agricultural subcatchment with texture-contrast soils.

Modelling
The ICECREAM model was first parameterised for the Nummenpää catchment.Next, a sensitivity analysis for the new macropore parameters was performed.Then, the model was calibrated for the reference period before the gypsum amendment through application of the monitoring data from the Nummenpää catchment.Next, the effect of gypsum amendment was added to the model in view of the various studies that were carried out during the TraP project and elsewhere.Finally, the model was verified against the monitoring data in the Nummenpää catchment for the period after the gypsum amendment.A detailed description of each phase (see Fig. 2) is given below.

Setting up the model
The model was set up to simulate average conditions for the three main agricultural land uses in the Nummenpää catchment: 1) spring barley with autumn ploughing, 2) spring barley with shallow cultivation, and 3) cabbage with autumn ploughing.Spring barley was chosen to represent spring cereals in general.The results of the three simulations were then combined on an area-weighted basis.The slope (1.6%) and soil type (silty clay loam) were kept constant in all simulations.Fertilisation rates were set at the average values used for barley (95 kg ha -1 N, 8 kg ha -1 P) and cabbage (225 kg ha -1 N, 70 kg ha -1 P) in the area.In addition, the average STP values for barley and cabbage fields were used.Simulations were run with 10-year climate data (2001-2010) from the Finnish Meteorological Institute, with the first year used for model spin-up.The date for seedbed preparation using a harrow, planting, and fertilisation was set to the beginning of May and that for harvesting to the end of August, approximating standard practice.Autumn ploughing (22 cm deep) and shallow cultivation (15 cm deep) was set to the end of October, as typical for the area.Both the efficiency parameters of tillage practices and the plant parameters of barley (see Appendix ) were adopted from earlier parameterisations (Bärlund andTattari 2001, Tattari et al. 2001).There are no existing modelling studies of P leaching from cabbage fields, so the parameters for cabbage were collected from different sources, such as Finnish agricultural fact sheets and the GLEAMS manual (Knisel 1993).
For cabbage, some of the parameters had to be based on expert judgement.
Model parameters related to soil characteristics were based on the article of Bärlund et al. (2009), in which typical soil profiles for mineral agricultural soils in Finland are specified.In the simulations, the field size was set to 2 ha (length 200 m, width 100 m), which is close to the Finnish mean field area, 2.2 ha (Puustinen et al. 1994) and is in accordance with earlier Finnish ICECREAM applications.The simulated profile was 1 m.The soil parameters are given separately for topsoil (0-25 cm; three layers) and subsoil (26-75 cm; four layers).

Sensitivity analysis of the macropore parameters
The parameterisation of the model runs for testing the sensitivity of the macropore parameters were based on the typical conditions in the Nummenpää catchment (e.g., spring barley, STP of 12 mg l −1 , slope 1.6%, silty clay loam, and autumn ploughing).
The model was first tested with macroporosity values of 5, 10, and 15% and with macropores always active when flow was generated (i.e. a threshold value of 0) and with the ω parameter (i.e. the percentage of PP from P pools entering macropores) set to 1•10 −5 .The amount of macropore flow and its share in the total flow increased linearly when macroporosity increased.The average proportion of macropore flow to total flow (i.e.subsurface + surface flow) ranged from 8% (5% macroporosity) to 24% (15% macroporosity).It turned out that even at 15% macroporosity, the matrix flow was greater than the macropore flow, with the average share of macropore flow in total subsurface flow being only 35%.Karvonen (2004) and Hintikka et al. (2008) have suggested that the percentage is as high as 87-99%.In ICECREAM, without the macropore flow, the annual subsurface flow was, on average, 61% of the total flow.However, when macropore flow was added to the model, the average percentage of subsurface flow out of total flow increased to 72%, which is in agreement with Turtola et al. (2007), who observed in a ploughed heavy clay field (slope 2%) that subsurface flow accounted for 58-92% of the total flow annually.The increase in the share of subsurface flow in ICECREAM is explained simply by the added macropore flow in which the water retention processes are lacking.The cumulative PP load in macropore flow was nearly constant across all three simulations, since the ω parameter does not depend on the amount of water flowing through the macropores (as long as there is flow).However, increasing the area of macropores decreases the share of surface flow in the total flow and thus affects P losses originating from surface runoff.With 15% macroporosity, a threshold value of 0, and ω being 1•10 −5 , the average PP concentration in the surface flow was 0.32 mg l -1 and in the subsurface flow 0.47 mg l -1 .The subsurface concentration was about twice the concentration found by Turtola and Paajanen (1995) for artificial drainage in a heavy clay field.However, in rainfall simulations performed by Uusitalo et al. (2012), the average concentrations of percolated PP in samples from ploughed fields were as high as 0.76 mg l -1 .Still, there are only a few studies (e.g., Uusitalo et al. 2001, Uusitalo et al. 2007) in which the concentration of different P forms has even occasionally been less in surface flow than in subsurface flow.
The model does not simulate any macropore flow if the amount of water reaching the soil surface is below 1 mm d -1 .In addition, macropores are activated only if the ratio of the actual water storage to the upper limit of the soil's storage is above the user-assigned threshold.The threshold value is soil-specific and depends on the water retention curve, which determines the upper limit for the threshold.With silty clay loam, this maximum threshold was 0.7, and the difference in cumulative PP in macropore flow between tested threshold values 0.0 and 0.7 averaged 0.63 kg ha −1 a −1 (see Fig. 3).
Calibration for the reference period The model was calibrated for the situation before the gypsum amendment in line with continuous monitoring data from 21 February to 21 May 2008 (i.e. the reference period).ICECREAM simulations were aggregated to cumulative curves for both DP and PP that included 15% cabbage (ploughing, STP 50 mg l −1 ) and 85% barley (59% of this from ploughing and 41% from cultivation, STP 16 mg l −1 ).The calibration was primarily based on the sensitivity analysis of the macropore parameters, since further efforts to calibrate the model with parameters potentially affecting results (e.g.snow melt parameters for the experimental area or the soil erodibility factor) did not improve the agreement between the simulations and the monitoring data.The area of macropores was set to 15%, as a higher percentage would exceed the probable percentage of macropore area for a clayey field (Aura 1995, Pietola and Yli-Halla 2007, Hintikka et al. 2008, Alakukku et al. 2010).With the threshold value 0.7 and with ω as 1•10 −6 , the model yielded the best simulation (see Fig. 5).Despite the calibration, ICECREAM was unable to predict the rapid increase in PP load during a rainy period starting on 5 April 2008 and lasting nine days, resulting in underestimation by 0.25 kg ha -1 after the three-month reference period for the simulated PP load.For DP, the calibration was more successful: After the three months, the simulated DP load was only 0.03 kg ha -1 higher than the load measured during monitoring.

Adding the effect of gypsum amendment
Since the ICECREAM model does not include chemical processes such as the dissolution of calcium sulphate, which increases ionic strength, other approaches had to be used to simulate the gypsum amendment.Although some experimental and laboratory data were available, the following assumptions were made on account of the scarcity of data: the effect of gypsum was assumed to be constant over time and independent of the slope and STP of the field.
In the laboratory experiments carried out for the soils from the Kisko farm, gypsum decreased DRP concentrations in percolation water.Over an STP range of 5-40 mg l −1 and without gypsum, DRP (mg l −1 ) in percolation water was related to STP in the relationship DRP = 0.0127•STP -0.0232 (r 2 = 0.85) (10) The gypsum treatment of 4000 kg ha −1 changed this relationship to DRP = 0.0084•STP -0.0583 (r 2 = 0.90) (11) For the STP range 16-50 mg l −1 (the average range in Nummenpää spring barley and cabbage fields), equations 10 and 11 would translate to a DRP decrease of nearly 50% due to gypsum application.The same effect was assumed to occur for DP, as typically a negligible proportion of the DP in Finnish clayey soils is present in a non-reactive form (Turtola 1996).In addition, Murphy and Stevens (2010) showed that gypsum can also decrease the solubility of organic P by 10-53%, which is often considered to represent the non-reactive DP.Accordingly, the equation linking STP and C Pil in the ICECREAM model was fixed such that C Pil fell to about half of the original concentration after the gypsum amendment and stayed constant thereafter.Fixing C Pil reduced DP, DP p , and DP bypass through equations 1-5 but also PP loss in surface water through the equation In addition, PP bypass decreased indirectly through Equation 12.This reduction in PP was welcome, since laboratory experiments by Uusitalo et al. (2012) indicate that gypsum affects PP loss even more than it does DRP leaching.According to their results, gypsum application resulted in lower PP and DRP concentration in percolating water, with this effect continuing after three winters from the gypsum application, with about 50% lower concentrations of TP as a result.monitoring seasons (I-IV), forced to start at the same level as the gypsum simulation so as to enable better According to Hämäläinen et al. (2010), soils treated with gypsum improved infiltration capacity under field conditions.This would suggest decreased surface runoff and thus lead to lower P loading on soils that are more prone to loading via surface processes.Also Cox et al. (2005) noticed that gypsum may enhance the share of interflow by increasing the size and stability of structural aggregates in the soil.One possibility for improving the infiltration in the model is to increase the saturated hydraulic conductivity.However, this increased percolation only negligibly.Therefore, the value of the SCS curve number was further lowered by five units, which resulted in decreased surface runoff and hence in lower mechanical surface erosion, which, in turn, increased the amount of macropore and percolation water.A decrease in erosion reduced PP loss in surface water through Equation 6 and PP bypass indirectly through Equation 12.In addition, the erodibility factor was lowered by 17%.
The lowering of the C Pil , soil erodibility factor and SCS curve numbers and increase in saturated hydraulic conductivity did not decrease PP loading enough, reducing the PP load by less than 40% and the DP load by around 50%.Therefore, further adjustment was made, by decreasing the amount of PP that entered macropores from the P pools by 20% (i.e., from 1•10 -5 to 8•10 -6 ).With this final change, the reduction of both PP and DP was close to 50%, which corresponded well with the prevailing laboratory studies.The changed parameters are shown in Table 2.

Model output and verification
The adjusted ICECREAM model was used to simulate the effect of gypsum amendment in Nummenpää for 5 November 2008 -18 June 2010 (see Fig. 6).The result was then compared with that obtained with the model in accordance with the parameterisation done for the reference period -i.e., the situation without the gypsum amendment.
The new (ICE0609) version simulated the effect of gypsum satisfactorily for three out of the four measurement seasons (seasons I-IV in Figure 6).According to the simulations, if no gypsum was applied, the cumulative TP load was 3.4 kg ha −1 for the entire period.With gypsum addition, the simulated load was 1.5 kg ha −1 (i.e.44%) lower.
Table 3 shows the simulated flows and P loads broken down into surface flow and subsurface flow.According to the model, before the gypsum amendment, 59-85% of water went through the soil profile.The changes made to the model to simulate the gypsum amendment increased the share of subsurface flow to 76-92%.Before the gypsum amendment, the simulated share of P load in surface flow out of total load was 17%.With gypsum, the figure dropped to 13%, meaning that gypsum not only reduced the P loss but also reduced the share of surface runoff P. According to our simulations, the gypsum amendment reduced both PP and DP nearly 50% during the high-flow seasons monitored.Comparison between the simulated and monitored TP load (see Table 3) shows that ICECREAM underestimates TP loading in all seasons except the third one, suggesting that the actual reduction by gypsum was slightly lower.If one ignores the third, poorly modelled season, both PP and DP varied at the same rate as the measured losses did.Over the entire modelled period, the simulated gypsum-induced reduction in TP load was 44% (see Fig. 6).Therefore, the modelled effect of gypsum on P losses was less during the nonmonitored, low-flow periods.

Discussion
Our main objective was to incorporate the effect of gypsum treatment into the ICECREAM model, which was here extended to include macropore flow.The effect of gypsum was added via rather simple modifications that were based on water analyses and other available information, since no detailed data for the soil processes after the gypsum treatment were available.The inclusion of macropore flow in the model made it a more reliable tool for quantifying P losses from fields.The ICE0608 version was not able to simulate PP in subsurface flow, as percolating water included only DP.In ICE0609, the share of PP in subsurface flow is better in line with studies (e.g.Uusitalo et al. 2001), which have shown that the dominant form of P in subsurface flow can be PP rather than DP.
The model calibration and validation was based on data from an outlet of the catchment that has diverse land use, while the model simulated only loading from fields.Not all sources of nonagricultural load could be accounted for, because of the limited data set for such factors as the actual amount, let alone temporal variation, of the sewage load from the scattered settlement in Nummenpää.As we were unable to account for it with a daily timestep, the calibrated model explains P loss from wastewater as well as from fields.Besides the varied land use, the monitoring data included the effect of in-stream processes (such as deposition, resuspension, and bank erosion).
Unfortunately such contingencies arising in the Nummenpää catchment were not studied in any way.However, the in-stream processes could explain the sudden increase in monitored PP load during the reference period, which the model was unable to simulate properly.Although different loading sources could have been taken into account with, for instance, catchment-scale models, addition of the effect of gypsum on agricultural soils would have been poorer on account of the more general field and soil descriptions.
Regardless of the hindrances mentioned above, the model was able to simulate gypsum addition in three of the four monitored seasons with reasonable accuracy.The overestimation of P loading in autumn 2009 may result from the fact that ICECREAM does not account for some processes that are central to the gypsum-induced pro-cesses, such as flocculation of the smallest soil particles.For now, lack of sufficient information on the flocculation processes prevents further changing and testing of the model in this respect.
Otherwise, the simulations yielded realistic results.The share of subsurface flow was close to the figures found in Finnish field experiments on ploughed fields, although the experiments also showed that the percentage can vary greatly, depending on soil, slope, cultivation method, the age of the drainage system, and weather (Vakkilainen et al. 2010).In addition, the changes that were made to ICECREAM to simulate the effect of gypsum amendment decreased the share of P load in surface flow out of total load by 4%.This is in agreement with the work of Cox et al. (2005), who noticed that gypsum appears to alter the relative proportions of P movement through overland flow and interflow.
Besides the gypsum-induced flocculation processes, there are other features to be improved in ICECREAM.Work on leaching through macropores could be developed further, for both spatial and temporal resolution.It now gives a constant amount of P (kg ha −1 d −1 ) through macropores, independent of the flow volume or the area of macropores.In addition, addressing of the relationship between matrix and macropores could be improved, and soil-specific time frames for initiation of macropore flow could be added to the model via further adjustment of the threshold value.In the Swedish ICECREAM model (Larsson et al. 2007), infiltration into macropores is also assumed to take place only if the water content of the upper layers exceeds a given soil moisture threshold.Authors involved in that work suggest that more physically based models could produce more accurate results.However, if the required parameter set is too difficult to measure, applying such a model becomes challenging.A literature survey showed that few data are available on nutrient concentrations in the water routed through macropores.These data would improve the parameterisation and render recoding of the macropore processes possible.To improve ICECREAM's applicability further, the calculation procedure for P sorption should be made better.Yli-Halla et al. (2005) noticed that ICECREAM does not fully cope with changes in STP at annual P application rates above 60 kg ha −1 , indicating that the immobilisation of large surpluses of fertiliser P in soil is stronger than currently predicted by the model.For such contexts as cabbage fields, application rates are usually higher.

Conclusions
The modified ICECREAM model enables comparison of performance between agrienvironmental measures in reduction of P fluxes under conditions similar to those of the Nummenpää fields.However, it is important to note that a model's capability of estimating the efficiency of a given water protection measure, especially a new measure, is rather unclear in the absence of sufficient data for proper calibration and testing of the model.Also, no verification of the model for other fields could be done in this study, because of the lack of data from other experiments wherein gypsum was spread and monitored, which would have improved the reliability of the modelling results.For accurate incorporation of the effect of gypsum into ICECREAM, or any P transport model, the flocculation and aggregation of the smallest soil constituents should be elucidated and the performance of gypsum in other than clayey soils should be studied in more detail.Additionally, more water quality data should be available, to allow proper estimation of the individual flow paths of nutrients.

Fig. 3 .
Fig. 3. Simulated cumulative loss of particulate phosphorus (PP) in macropore flow with threshold values of 0.0, 0.3, 0.5, and 0.7, indicating the soil moisture level in the two upper soil layers above which macropore flow starts.

Fig. 5 .Figure 5 .
Fig. 5. Calibrated model runs for the reference period (i.e., the period before the gypsum amendment) for cumulative particulate phosphorus (PP) and dissolved phosphorus (DP) losses, compared with the losses observed in the catchment experiment.

Fig. 6 .Figure 5 . 694 Figure 6 .
Fig.6.Simulated total phosphorus (TP) losses(5 Nov. 2008 -18 Jun.2010) plotted against the four monitoring seasons (I-IV), forced to start at the same level as the gypsum simulation so as to enable better comparison.'No gypsum' = cultivation practices assumed to stay the same; 'Gypsum' = amendment assumed to take place on 91% of the field area, as in the situation monitored in Nummenpää.

Table 1 :
Crops and cultivation methods in the Nummenpää catchment in 2009 Meteorological data (daily mean air temperature, precipitation, wind speed, relative humidity, and cloudiness) were collected from nearby stations (Helsinki-Vantaa Airport and the Nurmijärvi geophysical observatory) of the Finnish Meteorological Institute.Precipitation data were replaced with automatic weather station data for the periods for which they were available for the Nummenpää site(13-21 May 2008, 1 November 2008-14 January 2009,  13 March-16 June 2009, 31 August-18 December 2009, and 12 March-18 June 2010).

Table 3 :
Simulated and monitored flows and P loads (PP for particulate phosphorus and DP for dissolved phosphorus), the simulations being further divided into surface and subsurface flow -the division into seasons I -IV (see also Figure6) is due to periodic monitoring