Soil parameter variability affecting simulated field-scale water balance , erosion and phosphorus losses

Ilona Bärlund1,2, Sirkka Tattari2, Markku Puustinen2, Jari Koskiaho2, Markku Yli-Halla3 and Maximilian Posch4 1CESR, University of Kassel, Kurt-Wolters-Strasse 3, D-34125 Kassel, Germany, e-mail: baerlund@usf.uni-kassel.de 2Finnish Environment Institute, PO Box 140, FI-00251 Helsinki, Finland 3Department of Applied Chemistry and Microbiology, PO Box 27, FI-00014 University of Helsinki, Finland 4Coordination Centre for Effects,PBL, PO Box 303, NL-3720 AH Bilthoven, The Netherlands


Introduction
Erosion and nutrient losses are commonly modelled at the scale of a single field or at catchment scale (e.g.Rekolainen et al. 1999, Wade et al. 2004, Wolf et al. 2005, Grizzetti et al. 2005, Santhi et al. 2006, Bärlund et al. 2007, Hesse et al. 2008).Plot-scale erosion and nutrient loss data is often scarce and typically it is used to study the effects of different agricultural practices on environmental loading (e.g.Puustinen et al. 2005, Jégo et al. 2008).When longterm diffuse loading data is available, it is possible to get information on temporal variation of nutrient losses (Granlund et al. 2005, Puustinen et al. 2007).Field-scale modelling is an effective and useful approach to develop and assess the efficiency of management strategies which reflect natural processes, but the models are often poorly calibrated.When fixing the calibration of event-based modelling to a few experimental plots representing a small excerpt of possible soil-crop combinations we might lose variability that exists in nutrient losses at the catchment scale.The characteristics of fields, such as soil type, porosity, hydraulic conductivity, field capacity and wilting point, and the occurrence of horizons of different textures within a profile, can vary widely, even within a field and certainly within a catchment.The generalisation of modelling results based on a single field to a wide agricultural area may therefore result in erroneous conclusions about erosion and nutrient losses.
The parent material of Finnish soils was deposited quite recently (<12 000 years B.P.), and owing to rather weak soil formation in the cool climate, soils have been nationally classified according to texture and organic matter content, not according to pedogenesis.Soils are separated into (1) organic soils (>20% organic matter), (2) till soils (moraine, consisting of all or most textural fractions) and (3) sorted mineral soils, further divided into different textural classes.According to a compilation of soil test results of 700 000 topsoil samples and 37 000 subsoil samples collected in 1981-85 (Kähäri et al. 1987), silt prevails in the topsoils (share of soil classes 32%) while clay soils are most frequent in the subsoil (share 39%).Organic soils (18% of top-soils) are most common in the north and east; those with shallow organic topsoil have mineral subsoils of varying texture.Till soils (16% of topsoils) are common in east-central and northern Finland and clay soils on the southern and south-western coast while silt soils are more or less evenly distributed throughout the country (Kähäri et al. 1987).Commonly there is a clayey subsoil beneath a silty topsoil (Yli-Halla et al. 2000), the clay soils dominate in the major agricultural areas.Rivers draining these areas transport substantial amounts of eroded material into the sea, and the flow retention time between fields and coastal waters is short.
The ICECREAM model has been developed to study variables affecting erosion and processes governing phosphorus (P) and nitrogen (N) losses at field-scale in Finland (Tattari et al. 2001, Yli-Halla et al. 2005, Granlund et al. 2008).Results of single field-scale modelling have typically been extrapolated to describe effects of different management practices on agricultural land in a catchment (e.g.Mattila et al. 2007).This approach is based on modelling of different prevailing soil-crop-slope combinations in the study area.Geological soil maps, agricultural statistics and digital elevation data can be utilised to derive information on these factors but only on a very general level, usually without detailed information on single model input parameter values.
The most important question for model use remains: how to parameterise soil to catch some of the variability observed in soil sampling?A comprehensive study on physical properties of agricultural soils in Finland (KUTI, Puustinen et al. 1994) was carried out in the early 1990s.This data can serve to anchor the model parameterisation on a wider basis.In this study, typical soil profiles and the corresponding range of soil characteristics in Finland are described.Since no one correct parameter set exists, a sensitivity study was required to show the range of key parameters and the effect of this variability on simulated output variables.For this the ICECREAM model was run for a range of soil physical parameters.The output variables studied were hydrological variables as well as erosion and P loss.
The objectives of this study were (1) to parameterise typical mineral soil profiles for the ICECREAM model using key soil parameters, (2) to model water balance components, erosion and P losses for these typical Finnish agricultural soils to highlight those output variables that are especially affected by the variability in soil parameters and (3) to compare the results with measured data from field experiments.Moreover, the correlation of model results with model parameters was investigated.

Description of the soil profiles
This study utilised the data on physical properties of Finnish agricultural soils collected in the KUTI project (Puustinen et al. 1994).In that project the production infrastructure of the farms, the cultivation practices and farming methods in Finland were studied.One major aim of the KUTI project was to identify the field characteristics which affect erosion and nutrient losses, such as soil type, slope, plant available P and vicinity of the field to a watercourse or to a main ditch.Soil characteristics of the sampled fields (n=1065) were statistically analysed within the soil type groups used in Finland.The share of fields randomly sampled was weighted according to the share of agricultural land with this particular soil type.In each randomly selected study field, soil characteristics and other related variables were analysed from soil samples taken systematically along a transect across the study field.Results of the statistical analyses including averages, standard deviations and 5, 10, 25, 50, 75, 90 and 95% percentiles were calculated from all measured variables.According to the KUTI results, cultivated fields of Finland are generally quite flat, the median slope being 0.7-0.8%.The share of the field area with a slope greater than 3% is 17% (377 000 ha). Continuous authentic soil profiles were taken by auger sampler down to the depth of two meters and the soil type was specified visually in the field at 10 cm intervals and additionally in the laboratory from soil samples taken from the topsoil (plough layer, topmost 20-25 cm) and the subsoil beneath the plough layer.The share of clay soils in the subsoil of the KUTI material was 35% (770 000 ha), silty soils 16% (355 000 ha), coarse mineral soil 35% (775 000 ha) and organic soils 14% (317 000 ha).The KUTI material corresponds to the distribution of soil types in soil testing (Kähäri et al. 1987), with the exception that the coarse mineral soils were more abundant and silt soils less abundant in the KUTI material.
The soil profiles for the ICECREAM model were compiled from the KUTI database.The classification was based on the subsoil clay content.Six different mineral soil type groups (Groups A to F, three soils in each group) were composed on the basis of soil textures separately for the topsoil and the subsoil (Fig. 1).It is to be noted that the limits of silt (0.002-0.06 mm) and sand (0.06-2.0 mm) fractions as used in the U.S. classification, on which ICECREAM is based on, are different from the Finnish ones (0.002-0.02 mm and 0.02-0.2,respectively).Group A includes three clay soils (>30% of clay) representing (1) minimum, (2) average and (3) maximum clay contents within the clay soils (for further details see Table Ia, Annex 1).Groups B and C fall into the classes of silty clay loam, silty clay and clay in the subsoil, incorporating topsoils from silt loam to clay loam.Group D and E stand for silt and silt loam soils, group E being clearly more coarse-textured with an average sand content of 35% in the subsoil compared to only 5% for group D. The group F represents sandy loam and loamy sand and is clearly dominated by the sand fraction, all with clay contents <5% in the topsoil and <15% in the subsoil.The soil profiles constructed for the present study stand for typical conditions for Finnish cultivated soils where the topsoil is more coarse-textured than the subsoil (Mokma et al. 2000, Yli-Halla andMokma 2001).The three profiles of each group represent the range of soil textures in each soil group.Even though from Fig. 1 it seems that the profiles concentrate along the right side and the bottom of the classification triangle, it is estimated that these profiles cover approximately 70% of cultivated area in Finland.
Only the organic and the most coarse-grained soils are excluded in this study.
Porosity, field capacity, wilting point, saturated hydraulic conductivity in the subsoil and specific gravity are measured values of selected soils included in the KUTI database.Organic matter content was set as a constant for each soil group based on the measurements of selected soils in the topsoil layer and a unique value (1%) was given for all subsoils based on expert knowledge.Typical pH values (Lilja et al. 2006) were assigned to each soil group.Saturated hydraulic conductivity was only measured in the subsoil in each field studied in the KUTI project.The hydraulic conductivity of the corresponding topsoil layer was obtained by multiplying the subsoil value by four, based on selected data from MTT Agrifood Research Finland (Heinonen et al. 2001).The variability of topsoil saturated hydraulic conductivity over time is generally very large due to natural variability as well as tillage practices and machinery used.

The ICECREAM model
The ICECREAM model is based on the CREAMS model (Knisel 1980) but it has been further developed to better deal with the Finnish conditions (Rekolainen and Posch 1993, Tattari et al. 2001, Yli-Halla et al. 2005).The hydrology component simulates daily runoff using a modification of the SCS curve number method (Soil Conservation Service 1972).The erosion sub-model is based on the modified USLE (Foster et al. 1977).The P and N sub-model of ICECREAM is primarily based on the GLEAMS (Knisel 1993) and EPIC (Jones et al. 1984) models. Compared to Yli-Halla et al. (2005), the evapotranspiration is now computed with the Penman-Monteith model (Monteith and Unsworth 1995) instead of the original Ritchie equations (Ritchie 1972).
In ICECREAM the variability in soil physical properties affects directly water balance components, erosion as well as N and P loads.In addition, it affects parameters such as soil erodibility (K factor) in the USLE equation and the curve number CN2 used to calculate surface runoff.The K factor was calculated from soil texture, structure and permeability according to Knisel (1993).The CN2 value is an empirical parameter and the choice followed the outline used in GLEAMS (Knisel 1993, p. 47) classifying the different soil types according to the respective Hydrologic Soil Group for 'fallow' and 'small grain'.The Manning's n values are based on expert knowledge: for barley and stubble 0.040, for the use of harrow 0.022 and for the use of plough 0.025 were used for all soil types.The CN2 values and K factors vary with soil type and are thus presented in Tables Ib-c, Annex 1.
Previously (Yli-Halla et al. 2005), the initial state of the labile P pool (P il ) was operationally defined as P extracted by anion exchange resin (Uusitalo et al. 2001).An equation was developed to derive the P il from the results of soil test P, extracted with an acid (pH 4.65) ammonium acetate solution (P Ac ).The equation that was used in earlier versions of the model was based on analyses of 62 soil samples of different texture classes.For the present study, the equation was amended by analysing 91 more soil samples for the P extracted  1a in Annex 1).
with acid ammonium acetate (P Ac ) and anion exchange resin (P il ).These additional soils were predominantly clayey, representing the most important agricultural soils of southern and south-western Finland.The results of the two soil materials were combined, and the following equation ( 1) was obtained (R 2 =0.568):The range of P Ac values used in the calculation of equation ( 1) was 2.3-74 mg P dm -3 of soil.The relationship between P il and P Ac for P Ac ≤ 2.3 mg P dm -3 -a range that seldom occurs in field soils in practice -is expressed in equation ( 2), which is a straight line between 0 and the starting point of the P il value obtained from equation (1).The equation (1) has the same shape as the previous one (Yli-Halla et al. 2005), i.e. the higher the concentration of P il , the larger part of it is in the easily soluble form, extractable with ammonium acetate.The new equation yields slightly lower results for P il at any level of P Ac , a result compatible with the fact that at any level of plant-available P, fine-textured soils have lower soil test P levels.

Model set-up
All simulations were performed for a standard spring barley cultivation practice and 15-year climate data (1981-2000, with five years for 'warming-up' the model) from Jokioinen, south-western Finland.The standard cultivation practice means annual repetition of: harvesting at the end of August, ploughing early September and seedbed preparation using a harrow as well as sowing in early May.A fertilization of 15 kg ha -1 P and 100 kg ha -1 N was added each year at sowing.A field slope of 3%, which is higher than the average field slope, was selected because differences in simulated results are better visible, and 3% is not an extreme value.In the simulations, the field size was set to 2.2 ha (length of the field 183 m, width 120 m), the simulated profile thickness was 1 m (topsoil 22 cm, subsoil 78 cm) and the initial soil test P Ac to 12 mg dm -3 , both being average values in Finland.The model was run for each soil profile (see Tables Ia-c, Annex 1) separately and the simulation results presented are averages of annual sums over the 15-year simulation period.

Description of the statistical method
The correlation matrix between simulated erosion as well as phosphorus losses and the initial model soil parameters were calculated with the SPSS® software.Here, the analysis was first made for all six soil type groups and in addition the soils A, B and C (clayey soils) and soils D and E (silty soils) were analysed separately.Only parameters that were measured in the KUTI project were included in the analysis.Thus soil organic matter content, pH, CN2, the erodibility factor (KSOIL) and Manning's n were excluded.The statistical analysis was first and foremost performed to study how the model parameters correlate with model output in general.Since the correlation coefficient r is a rather poor statistic measure for deciding whether the observed correlation is statistically significant, no further conclusions can be drawn for the causal relations of the model parameters and the output.

Results and discussion
The variability in the output depended on the soil type and differed from one output variable to another.One reason is the variability in the measured input parameters describing soil texture, such as the difference in soil moisture content at field capacity and at wilting point (i.e.available water capacity), where the average is largest for the medium-textured soils of profiles D and E but with a considerable variability also for especially the topsoil of profile C (Fig. 2).Also the variability in saturated hydraulic conductivity was remarkable, both between the profile groups and within one group, especially for the coarsest profiles E and F. From previous sensitivity studies (Bärlund and Tattari 2001) it seems, however, that a model with a water-holding-capacity-type strategy to store and transport water is much more sensitive to the values of the pF curve than to the saturated hydraulic conductivity.
The average annual precipitation in 1986-2000 was 687 mm a -1 .On an annual level 40-55% of precipitation was simulated to be either transpirated from vegetation or evaporated from the soil depending on the soil profile.The highest average evapotranspiration was for soil profile E (363 mm a -1 ) and lowest for soil profile C (278 mm a -1 ).The total simulated evapotranspiration values are of the same magnitude as or slightly smaller than presented by Solantie and Joukola (2001) for agricultural fields in southern Finland.Simulated evaporation was lower than transpiration for all soil profiles (Fig. 3).The variation of transpiration between the six soil profiles was considerable (max.68 mm a -1 ), but the variation within a soil group was relatively small.On the other hand, the variation of soil evaporation between the soil profiles was small (22 mm a -1 ), but here the variation within the groups was rather large, being highest for the medium soil groups C-E where also the variability in the parameters of the pF curve was highest.One of the main factors affecting transpiration and evaporation was the available water capacity of soil.Evaporation was stronger linked to the actual soil water content whereas the uniform crop reduced variability in transpiration.One reason for the highest variability in evaporation in the soil groups C-E was linked to the highest variability in CN2 for these groups which affected the amount of infiltrated water and thus the actual soil water content.
The relationship between surface runoff and root zone percolation showed a distinct pattern of capacity type soil water management: the lower the surface runoff the higher the percolation (Fig. 4).The soil profile group E showed, however, very similar values to group A in spite of very different clay contents and pF curve values in these two groups.Soil type E has the highest available water capacity combined with the highest evapotranspiration values.This is one reason for the relatively low percolation compared to surface runoff, as percolation is the last calculated variable in a time step reflecting changes in all other output variables.The variability was large within all soil type groups but especially for groups A, B, D and E. These are groups with high variation in the CN2 values and the variability in surface runoff translates directly to the variability in percolation.This result highlights the importance of the empirical CN2 parameter.It can be adjusted to correspond with the actual drainage efficiency -the poorer the drainage, the more surface runoff is created and the higher is the CN2 value to be set.The ICECREAM version used did not include macropore flow.The effect of drainage improvement on the division between surface runoff and subsurface drainage has been demonstrated in the clayey Kotkanoja research field in Jokioinen, Finland by Knisel and Turtola (2000).As the aim of this study was to establish a general soil physical parameter set that can be used as a basis for simulations in whole of Finland, this site specific parameter was adjusted in a conservative way assuming average drainage efficiency.Simulated erosion followed to a certain degree the result of surface runoff but also the combined effect of soil texture factors (silt content) and factors derived from them (K factor) played a role in explaining the variability within the soil type.The annual average erosion values ranged from 32 to 563 kg ha -1 a -1 (Fig. 5).In an experimental field similar to soil A the measured annual average (1991-2001) erosion (surface plus drainage) for an autumn tilled field with changing drainage efficiency was ca.890 kg ha -1 a -1 (Turtola et al. 2007) which is somewhat higher than the range of variability of the simulated erosion.The simulated erosion, however, only accounts for erosion and sediment loss through surface runoff.The most erodible soil types were those of profile groups D and E. Both groups also showed by far the highest variability.
The results on P losses in surface runoff followed the results of surface runoff and erosion to a certain degree (Fig. 6).Except for the soil groups E and F, the average values of soluble P in surface runoff were very similar, around 0.27 kg ha -1 a -1 .In soil group E the very low clay content in the surface layer is conducive to a low P sorption capacity.This soil feature, combined with somewhat higher surface runoff, leads to higher dissolved P losses at the surface.Clearly the most important driving factor for particulate P loss in surface runoff was the simulated erosion, displayed also in the variability within the soil groups.Transport of particulate P through the drainage is still a transport pathway that is missing in the current ICECREAM version.This state of the art reduces the possibilities to compare the simulated results with field measurements on clay soils where losses of particulate P through subsurface drainage are considerable (Turtola and Paajanen 1995).An approach to include macropore transport in fieldscale simulation of P loss is given e.g. by Larsson et al. (2007).The soluble P loss in root zone percolation (Fig. 7) was closely connected to the hydrological variable, i.e. deep percolation (Fig. 4).Compared to the surface runoff component the amounts of P transported out of the profile were very small.The variability in soluble P loss in percolation reflected the variability of e.g. the clay content affecting sorption capacity but the differences between the soil groups were governed by the deep percolation result which became evident for the soil group E. Simulated P uptake of the total crop was higher than P fertilization (15 kg ha -1 a -1 ) for all soil groups (Fig. 8).The variability between the soil groups was small ranging from 20.5 to 21.3 kg ha -1 a -1 but the variability within a soil group (e.g.soil group C) can be as large.The main factor affecting the variability between the soil groups was transpiration (Fig. 3) but most likely also the amount of plant available P in the zone played a role as it is affected by the sorption capacity.Soil P pools were however not included in this analysis.The measured total P losses for a field similar to soil B (Turtola, 1999) were of the same order than the ones simulated in this study.The measured dissolved P loss was lower than the simulated one due to the higher P status in the model set-up.
The correlation matrix between simulated erosion and P losses and soil parameters is presented in Tables 1-3.Special emphasis is given here on r values > 0.5.For erosion the closest correlation was found with the fraction of sand and the saturated hydraulic conductivity.Additionally, for the silt soil group (Table 3) the correlation between erosion and the fraction of clay and sand as well as with available water in the plough layer was rather high.Considering all soil type groups together, the loss of particulate P followed the behaviour of erosion (Table 1).For the clay soil group (Table 2), saturated hydraulic conductivity was the only soil characteristic correlating well with particulate P. On the other hand, for the silt soil group several parameters, such as the fraction of clay, correlated positively and the fraction of sand and saturated hydraulic conductivity correlated negatively with particulate P loss.None of the parameters correlated very closely with soluble P in surface runoff (DPsurf), even though for clay and silt soil group the correlation was slightly better.In reality, soluble P transport is regulated by the soil P status and the cultivation practices.On the contrary, most parameters correlated quite closely with soluble P in root zone percolation (DPperc).Although the correlation was high, the amount of soluble P transported by this route, i.e. through percolation, was negligible (see Fig. 7).

Conclusions and outlook
This study presents parameter sets of soil physical properties for typical Finnish agricultural mineral soils that are based on existing soil data and their systematic use for model parameterisation purposes.
The aim was to generate a parameter set that can be used for field-scale modelling when the target is not to depict one particular field but e.g. a Finnish silty clay soil in general.The comparison of ICECREAM results with experimental data on the water balance and P losses is not to be considered a model validation, but to show that with the parameter sets proposed here a realistic magnitude of these selected output variables can be obtained.Further work with respect to nitrogen leaching is still required.
In this study two sources of variability were detected.The principal variability has its origin in the multitude of texture combinations that make one particular soil type.The second, sometimes even larger variability arises from the method to derive other parameters than the measured ones, such as the erodibility K factor and the CN2 values which are partially derived from soil physical parameters.On a field slope of 3% choices affecting erosion and thus particulate P transport in surface runoff are particularly important.
The result of the correlation analysis confirms the well-known close connection between erosion and particulate P transport.There seems to be a distinct difference, however, between the two soluble P loss variables.Soluble P in percolation was closely connected to the soil texture parameters included in this database, whereas for soluble P in runoff other parameters such as soil P status and the cultivation practices could be more relevant.Concerning soluble P loss, the silt loam soils (soil group E) posed an exception.As the soluble P loss in surface runoff for the other soil groups was nearly the same (and thus most likely not affected by variability in soil texture), for soil E texture seemed to play a role.For soluble P loss in percolation for group E the result was more affected by the indirect effect through reduced percolated water (resulting from enhanced evapotranspiration) than the effect of clay content on the sorption capacity.This highlights the complex, often non-linear relationships between the different processes and output components in a model like ICEREAM.
This study shows effects of parameter variability on simulated P transport.To cope with the ef- In this study, the basic parameterisation for the soils in Finland was established for modelling purposes.In the next phase, different classes of slope steepness will be added to simulate erosion and P losses using the model parameterisation presented here.The current ICECREAM work to create a parameter database for general model applications was limited to the physical parameters of soil.The parameterisation of crops would be the next step recommended for systematising the use of ICECREAM and similar type of models in Finland.

Fig. 1 .
Fig. 1.Clay, silt and sand fractions for the six soil type classes (A-F).Mean (clay) values are indicated as full circles for the top soil and open circles for the subsoils (connected by a dashed line).The arrows indicate the ranges used for each class in the simulations (A1 and A3, etc.) connected by thin dashed lines (see Table1ain Annex 1).

Fig. 2 .Fig. 3 .
Fig.2.Available water capacity (soil moisture content at field capacity minus soil moisture content at wilting point, average with min-max range) for the topsoils and subsoils of the typical soil profiles.

Table 1 .
Correlation of the KUTI-ICECREAM parameters with erosion phosphorus fractions for all soil textures.

Table 2 .
Correlation of the KUTI-ICECREAM parameters with erosion and phosphorus fractions for clay soils (soil groups A−C).

Table 3 .
Correlation of the KUTI-ICECREAM parameters with erosion for silt soils (soil groups D and E).

Table Ib .
Parameter values for the various soil type classes (A−F) and the minimum, mean and maximum clay content profiles within them (1−3).Parameters 7−12.

Table Ic .
Parameter values for the various soil type classes (A−F) and the minimum, mean and maximum clay content profiles within them (1−3).Parameters 13−18.Organic matter content for subsoil (SOLORG_j) of 0.01 m 3 m -3 was used for all soil types.