Appendix 2. Model algorithms and parameters used in PSYCHIC
Eunice Lord1, Mark Betson1, Paul Davison1, Phillip Owens2, Gavin Wood2 and Des Walling3
1ADAS Environment Systems Group, ADAS Wolverhampton.
2National Soil Resources Institute, Cranfield.
3Geography Department, University of Exeter.
1.Introduction
This Appendix contains a technical description of the PSYCHIC model: the algorithms used within the model, the values of relevant parameters and a user guide to running and operating the model. As previously noted, PSYCHIC is a process-based model, operating at both field and catchment scale. Although the underpinning science is the same in both cases, there are a number of assumptions built into the catchment scale model to allow necessary parameters to be inferred from available national datasets. The field scale model contains no such assumptions and therefore requires more input from the user than the catchment scale model. Below, we describe the operation of the catchment scale model. In most cases where the field scale model differs in input requirements, the difference is simply that the user is required to supply directly the parameter which, at catchment scale, is inferred.
2.Model algorithms
2.1Soil properties
In order to accurately estimate the relative importance of different hydrological pathways in a given field or grid square, it is essential to know a number of properties of the dominant soil. Unfortunately, such data are often not available, and it is necessary to infer them from such parameters as are available. PSYCHIC infers soil properties from soil series and land use.
The look-up-table (LUT) from which most soil properties are inferred, based on soil series, is given in Appendix A. It can be seen that one of the inferred properties is soil HOST class. This is then used, in conjunction with land use (grass or arable) to infer the presence (or not) of tile or mole drains (Appendix B). Other soil properties presented in Appendix A include soil particle size distribution (sand, silt and clay), organic carbon (OC) content, bulk density (BD), soil cohesivity (COH), soil total P (TP) and Olsen P. Organic carbon content and bulk density are functions of land use, with individual values for arable land, grass and “other” (e.g. woodland). It will be seen that there are some missing data in Appendix A: where the required value of bulk density or organic carbon content is missing, the “arable” value for the appropriate soil series is used.
2.2Water Balance
Within the PSYCHIC model is included a version of MCDM (Anthony, 2003) that calculates a water balance for seven representative land use types, bare soil, winter wheat, spring barley, potatoes, grass, rough grass and forest. MCDM is a modelling framework intended to calculate monthly values of potential and actual evapotranspiration, soil moisture deficit and soil drainage. It uses as input actual mean monthly climate values or if these are not available a default set of monthly values generated from long-term mean climate data (CRU 1961-1990) together with the choice of land use type and soil textures for a UK location defined by easting, northing and altitude. The model also predicts mean soil temperatures at a reference depth and observation hour using the homogenous conductor model (de Vries 1963, Gilman 1977).
2.2.1 Evapotranspiration
The method of calculating potential evapotranspiration is based on the Penman-Monteith equation (Monteith, 1981), using UK specific parameters taken from the MORECS and IRRIGUIDE models (Hough and Jones, 1997, Bailey and Spackman, 1996). Actual evapotranspiration is calculated using a modified version of the Thomas monthly book-keeping methodology (Thomas, 1981).Modifications have been made to allow for efficient evaporation of water intercepted by the crop canopy, prior to evapotranspiration from the soil moisture store (Alley, 1984).
2.2.2Surface Runoff
Surface runoff within PSYCHIC is calculated based on simple functions derived through meta-modelling the process driven Infiltration Excess model (IE, Anthony 2003). The IE model is a system for estimating the percentage of rainfall that becomes infiltration excess runoff, as a function of rainfall intensity and soil hydraulic conductivity.
The distribution of rainfall volume on rain days is modelled using a gamma distribution with the same regression analysis as Davison et al. (2005). The distribution of rainfall intensities for a daily rainfall volume is modelled using an exponential distribution of the cumulative depth of rainfall falling at a given intensity, sensu van Dijk (2002). This assumption produces results comparable to the intensity-duration studies of Bilham (1936). The parameters of the exponential distribution are estimated by an iterative process that compares the integrated kinetic energy of the storm (Marshall and Palmer, 1948) with empirical estimates of the mean expected energy for the storm volume, derived from regression analyses of 10 years of 15 minute rainfall data at 11 sites in England and Wales (Hutchins, 2002).
Storage of excess rainfall is calculated using the empirical method of Onstad (1984), as a function of surface random roughness and slope. This is the method used by WEPP (Flanagan and Nearing, 1995). Random roughness can be estimated from the roughness ratio coefficients used by EUROSEM, using an equation developed by Auerswald (1992). Surface runoff is calculated using a spatially variable infiltration model (Yu, 1999), based on the assumption that an exponential distribution of maximum infiltration rates occurs in the landscape. At increasing rainfall intensities, an increasingly large proportion of the landscape starts to contribute to runoff, as expressed by an exponential relationship between rainfall intensity and the area mean infiltration rate at that intensity.
Infiltration is calculated using the Green and Ampt (1911) model, and a mean calculated across the distribution of infiltration rates. The soil parameters required by the model are estimated from the particle size distribution, the organic carbon content and the bulk density using the HYPRES pedo-transfer functions (Nemes et al., 1999), and application of the Mualem van Genuchten equation (van Genuchten et al., 1991).
2.2.3Meta-function derivation
IE produces estimates of the percentage of monthly rainfall that becomes runoff. The first step in meta-modelling this process was to calculate the percentage runoff generated by a numerical integration of the two parameter gamma rainfall distribution above a given infiltration acceptance value. We fit a simplified function to predict the percentage runoff given the acceptance value (equation 1, Fig. 1).
1.
Where Acc is the acceptance capacity (mm), M is the rainfall per rain day (mm) and RO is the percentage rainfall runoff.
Figure 1. Acceptance parameter (mm) as a function of percentage runoff.
Using the assumption that the soil for each month starts at field capacity the IE model was run for 351 soil series using the input parameters of percentage sand, silt and clay, bulk density and the percentage organic carbon for arable soils. Representative rainfall and number of rain days were selected for four climatically different locations in England and Wales (Table 1) using data from the CRU 1961-1990 database (Barrow et al., 1993).
Table 1. Rainfall amount and days for selected climatic regions.
East / North / Location Description / Total annual rainfall (mm) / Total annual raindays (n)600000 / 290000 / East Anglia / 626.3 / 161.5
240000 / 800000 / Devon / 1227.0 / 200.2
270000 / 340000 / North West Wales / 2393.4 / 236.5
450000 / 140000 / Wiltshire / 849.3 / 166.8
The IE model estimate of percentage runoff per month for each series was converted to an acceptance parameter using equation 1. Comparing M to the acceptance parameter indicated a strong linear relationship for each soil series (Fig. 2).
Figure 2. Acceptance parameter (mm) as a function of rain per rain day (mm)
Using calculated values of M for all of the England and Wales locations a linear regression was performed on each of the soil series to produce values for the gradient and intercept. Finally the gradient and intercept values were compared against the original percentage sand, silt and clay values for each of the soil series to indicate any direct relationship. The result of this indicated a good linear relation to the silt content (Z) of the soil for both the gradient and the intercept (equations 2 and 3, Fig. 3).
2.
3.
There is also a strong correlation between intercept and the sand content of the soils (Fig. 3). An assessment of the correlation between the five soil parameters that drive the IE model and the linear regression gradient and intercept indicated the relative significance of each as a predictor for the acceptance value (Table 2).
Table 2. Correlation coefficients for the five key soil parameters for IE against the gradient and the intercept of the linear regression of the acceptance value and percentage runoff.
Sand / Silt / Clay / Bulk Density / Organic CarbonGradient / 0.850 / -0.949 / -0.321 / 0.122 / -0.021
Intercept / 0.838 / -0.950 / -0.297 / 0.113 / -0.014
Figure 3. Relationship between silt content and the gradient and intercept parameters obtained from the Infiltration Excess model.
As was suggested by Figure 3 silt is the strongest predictor followed by sand. Although silt is a strong enough predictor to yield an r2 value of 0.9 for both the gradient and the intercept, multiple regression analysis was undertaken by combining the silt parameter with the next strongest predictor to try to improve this (Table 3).
Table 3. The r2 values for the multiple regressions of the five soil parameters against the regression coefficients. *Inclusion of the clay parameter invalidated the regression tool being used.
Test / 1 / 2 / 3 / 4 / 5Silt (Z) / Yes / Yes / Yes / Yes / Yes
Sand (S) / No / Yes / Yes / Yes / Yes
Clay (C) / No / No / Yes / No / No
Bulk Density (BD) / No / No / No / Yes / Yes
Organic C (OC) / No / No / No / No / Yes
r2 (Grad.) / 0.90 / 0.90 / * / 0.94 / 0.98
r2 (Int.) / 0.90 / 0.90 / * / 0.94 / 0.98
The multiple regression tool within the Statistica software package was used for the analysis and within this it was found that the incorporation of the clay parameter invalidated the limits of the tool and therefore it had to be omitted from the results. The results of the multiple regression indicate that the most improvement from single parameter silt in the explained variance can be achieved using the silt, sand, bulk density and organic carbon; use of fewer parameters does not significantly improve the r2 value. The revised regression equations (4 and 5) are therefore,
(Equation 4)
(Equation 5)
2.2.4 Soil compaction, organic carbon content and water deficit
Equations 4 and 5 have been calculated using IE values determined for arable soils initially at field capacity. However soil moisture deficits can develop in soils during summer months and the OC contents of grassland soils are different from those of arable soils. The HYPRES equations within IE are sensitive to soil OC content and therefore will produce different results for grassland soils with a higher OC value. The NSRI soil series database incorporates four land use types; arable, ley grassland, permanent grassland and other land use with characteristic OC values for each. The analysis above has been repeated for the 314 soil series in the permanent grassland database as a representative OC value for non-arable land use types.
The resulting four parameter functions for the gradient and the intercept have been compared to the original functions (equations 4 and 5) for five characteristic soil series to indicate if there was a significant difference in the predicted acceptance value (Table 4).
Table 4. Comparison of the acceptance (Acc) parameter for five characteristic soil series predicted with the relations derived from arable and grassland soil types.
Series / Acc (arable) / Acc (grass) / % DifferenceCUCKNEY / 6.80 / 6.75 / 0.71
DENCHWORTH / 4.85 / 4.94 / 1.81
HANSLOPE / 5.25 / 5.33 / 1.53
MIDDLETON / 4.21 / 4.20 / 0.16
WICK / 5.97 / 5.95 / 0.41
In addition to the OC differences related to land use significant runoff differences are expected for compacted soils. To represent this the bulk density of each of the soil series has been increased by 35% based on work published by O’Sullivan et al. (1999) and experimenting with the simplified model of soil compaction they present.
The gradient and the intercept for both arable and grassland soils were recalculated with the increased values for bulk density and the corresponding multiple regression equations defined. Again the results using these relations with the five test soil series were compared with predictions using equations 4 and 5 with the increased bulk densities to assess if there was any significant improvement (Table 5).
Table 5. Comparison of the acceptance parameter for five characteristic soil series predicted with the relations derived from arable, compacted arable and compacted grassland soil types.
Series / Acc (arable) / Acc (compact arable) / %Diff. / Acc (arable) / Acc (compact grass) / %
Diff.
CUCKNEY / 3.87 / 3.44 / 11.21 / 5.13 / 4.84 / 5.69
DENCHWORTH / 3.00 / 2.74 / 8.85 / 3.50 / 3.29 / 6.00
HANSLOPE / 3.48 / 3.15 / 9.42 / 3.71 / 3.47 / 6.46
MIDDLETON / 1.97 / 1.96 / 0.46 / 2.55 / 2.45 / 4.17
WICK / 3.49 / 3.19 / 8.58 / 4.40 / 4.16 / 5.37
The comparison in Table 4 indicates no significant difference between the acceptance values calculated via regression analysis of IE results for grassland soil types and that from the existing functions derived from arable soil types (<2%). A similar approach for compacted grass and arable soils compared to the results using equations 4 and 5 produced more significant differences (<12% for compacted arable and <7% for compacted grass), Table 5.
Based on these results it was recommended that the derived meta-functions (equations 4 and 5) be used with modified soil parameters for grass and compacted soils rather than switching to specific meta-functions derived for these conditions from IE. In the case of grassland soils, the difference produced is not significant and for the compacted soils uncertainty about the response of IE to artificially high BD values for compacted soils means that it is difficult to justify using the specific functions produced.
The acceptance parameter for soil infiltration is expected to significantly increase with the development of a soil moisture deficit (SMD). A function has been derived to modify the calculated acceptance parameter in the presence of a SMD, such that acceptance reaches a maximum value (= (+1) ACC) ) during periods of high SMD and tends to the originally calculated value as SMD tends to zero. The function varies exponentially between the two extremes (equation 6). Another parameter () controls the shape of the exponential curve; if is low there is a rapid response to an increase in SMD and if high a reduced response.
6.
Note the value of SMD cannot equal zero in equation 6, it is assumed that if SMD = 0, ACCSMD = ACC. The values for and will depend on expert judgement in the absence of observed or published results, but as a first approximation values of = 1 and = 14 appear to significantly attenuate runoff where any SMD occurs (Figs. 4 and 5).
Figure 4. Predicted runoff for a site in North West Wales (see Table 1).
Figure 5. Predicted runoff for a site in East Anglia (see Table 1).
If the gradient defined by equation 4 is m and the intercept by equation 5 is c then combining these with equations 1 and 6 gives the following relation describing the percentage runoff.
7.
The calculated runoff (mm) in a given month is thus:
,8.
where Qm is monthly total runoff (mm) and RM monthly total rainfall (mm).
2.2.5Flow to tile or mole drains
As already discussed, MCDM diagnoses the water flow along each major hydrological pathway for each of the PSYCHIC land use classes for each month of the year. Included in these pathways is subsurface lateral flow, which notionally includes shallow lateral flow through the soil matrix and, if the landscape is underdrained, flow through drains. The method used to infer the presence of drains has been described in section 2.1.
A fraction of the diagnosed total subsurface flow is assumed to pass through tile or mole drains based on the HOST class of the drained soil. HER is assumed to pass along various pathways (overland flow, vertical flow, bypass flow, shallow lateral flow and deep lateral flow) in the proportions given by Boorman et al. (1995). The fraction of HER going to drainflow is determined as the ratio of shallow lateral flow to total vertical flow.
2.3Applications of P to agricultural land
Rates of application of manure (and hence phosphate) to agricultural land are derived from the Manures Management Database (Comber et al., 2005), to which the reader is referred for a full description of the database and the algorithms therein. A brief description of the major points of the database follows.
The database uses agricultural census data to calculate the area of managed grassland, rough grassland and arable land in each 1 km2 in England and Wales. The same source data are used to calculate the number of animals in each grid square in each census category. Export coefficients are applied to calculate manure production by animal type for each grid square, and look-up-tables (LUTs) used to determine the fraction of total manure of each type (eg. FYM, slurry) and the available (soluble) P in each manure type. The apportionment of manure between different types varies between counties.
Occasionally, “hotspots” in manure production are observed. It is thought that these are probably due to anomalous census returns, for example a farmer returning his herd as being entirely resident in the parish where he lives, when in fact they may be spread over a considerable surrounding area. Any such hotspots are removed by conservatively distributing the animals among surrounding grid squares.
Total manure production is then divided into applications to arable land and grass. Voided manure (i.e. that from grazing animals) is obviously applied only to grassland, while managed manure is applied to arable and grass in proportion to the area of each in the grid square. Finally, the number of manure categories is reduced by aggregating similar manure types to produce the following datasets, each of which contains the quantity of soluble phosphate applied to arable and grassland for each month of the year, in each grid square: