Inverse method to estimate mineralisation rate constants for nitrogen simulation models: interaction between sampling strategy and quality of parameter estimates.
Excessive nitrate leaching is an ever-increasing concern in Australia and elsewhere. This concern arises from both the contribution of nitrate leaching to soil acidification (Helyar and Porter 1989), and the contamination of groundwater (de Willigen and Ehlert 1996). Many land management practices lead to rates of nitrate leaching that are considerably greater than those occurring beneath `natural' ecosystems. These land management practices include the application of wastes from livestock production onto pastures or crops (Popp and Wichmann 1996; Chapman 1996; Cook et al. 1996). There is increasing pressure for land application of wastes (ARMCANZ and ANZECC 1995), many of which have high nitrogen contents and a high potential for nitrate generation and leaching. Sustainable land application of wastes requires that the accession of nitrate to groundwater be within acceptable limits. Development and testing of management practices that achieve this goal is most effectively done by targeted experimentation combined with modelling the fate of nitrogen in the wastes.
Prediction of the fate of nitrogen is accomplished by the use of a nitrogen cycling model coupled with models to simulate the response of plant growth and water movement to climatic variables. The nitrogen cycle is complex and demands a complex model with several rate parameters to simulate adequately the fate of nitrate. Successful application of these models requires unbiased estimates of the parameters describing the various processes within the nitrogen cycle.
Obtaining estimates of these model parameters can be an onerous task. Through judicious experimental design many parameters can be measured or estimated independently of the model, or they can be assumed of suficciently small importance that they can be assigned a default value. Some parameters, however, will normally be left undetermined. One technique for estimating these parameters is by using inverse methods based on field sampling and model optimisation.
Although the net nitrogen mineralisation flux can be determined by one of several methods, including various laboratory and field incubation methods (see Jarvis et al. 1996; Raison et al. 1987), difficulties arise in the estimation of the model's mineralisation rate parameter from mineralisation flux data. Most simulation models use first-order mineralisation submodels and partition the organic pool into several fractions that have different turnover rates (e.g. Johnsson et al. 1987; Metherell et al. 1993). Consequently, the mineralisation rate constant cannot directly be inferred from mineralisation flux data as it is convoluted with the partitioning of soil organic matter between the various pools.
Both the partitioning and the rate parameters are model-specific so, in general, either the model, or an analytical version of the mineralisation submodel, should be applied to the data in order to estimate the parameters. When applying the model to field conditions, the partitioning of the total soil organic matter is supplied to the model as an initial condition. With simulation time, organic matter distribution between the pools adjusts and becomes insensitive to the initial values supplied. Because incubation experiments are generally run for comparatively short times, the initial partitioning of the soil organic matter has a large effect on the model prediction of mineralisation flux.
Furthermore, mineralisation measurement techniques that generate conditions different from those normally found in the field incur a risk of unacceptably changing the parameters that are to be estimated from the data. For these reasons there is considerable merit in performing a calibration under field conditions using optimisation techniques employing the simulation model (Yamaguchi et al. 1992; Larocque and Banton 1996).
Inverse, or calibration, experiments used to estimate the nitrogen mineralisation parameters normally monitor nitrate concentrations in the soil or soil solution, at specific times and depths, in response to application of a source of nitrogen (e.g. waste or inorganic fertiliser). The water regime is usually also monitored. Following a spike addition of the nitrogen source, nitrate concentration at a given depth normally would be low initially, rise as the pulse of nitrate passed that depth, and then fall to background levels. Where the nitrogen source must first mineralise to produce nitrate, the appearance of the first rise in concentration may be delayed. Additionally there may be several peaks resulting from the effects of changes in the soil moisture and temperature on nitrogen mineralisation.
Logistics often preclude intensive sampling, and sampling may even cease before the passage of the nitrate pulse is complete. This results in uncertainty in the shape of the nitrate pulse and error in the calculation of the mass of nitrate leached. Added to this error is the inability of some sampling devices to access all of the soil water and solute (Brandi-Dohrn et al. 1996). Spatial variability and analytical errors add further uncertainty.
This combined error in the measurements can lead to bias and poor precision in the parameters estimated by inverse methods. The rate parameters are usually obtained from experiments carried out over a comparatively short time interval but will be used for simulations over a much longer time. It is therefore important that the parameters are unbiased and precise; it is particularly crucial that the mass of nutrient leached be accurately simulated.
Under ideal conditions, it may be possible to obtain samples of the composition and amount of leachate at intervals close enough in time to determine accurately the pattern and quantity of leaching to groundwater. If the experiment is carried out under rain-fed conditions, the time course of drainage and nitrate leaching may be highly episodic. This will contribute to the difficulty of adequate resolution of nitrate concentration sampling and will also adversely affect the uncertainty of water flow predictions by the model. Both of these compromise the success of the inverse technique for estimating parameter values.
Many of these concerns could be addressed adequately if a weighing lysimeter with automated sampling is used to monitor leachate quantity and quality. These devices are, however, expensive and hence are rarely available. Furthermore, they can be difficult to maintain, and there is a reluctance to destructively sample the soil inside them. Recently, simpler devices which provide the ability to measure the mass flux of a solute have become more routinely available. Such devices include drainage lysimeters (e.g. Cameron et al. 1992; Moyer et al. 1996), and mini-lysimeters with ion exchange resin to trap the solute (DiStefano and Gholz 1986; Schnabel et al. 1993; Sakadevan et al. 1994). The cost of their simplicity is usually that there is less temporal resolution of the pattern of the leaching.
This study investigates the potential benefits of the data collected using the simpler lysimeter designs for reducing bias and improving precision of estimates of nitrogen mineralisation parameters. Using numerical simulation, we evaluate the effects of sampling frequency and spatial variability in nitrate concentrations on parameter estimation, and demonstrate the improvements in parameter precision and unbias that result from adding information about the mass of nitrate leached and the water regime.
Leaching of nitrate following the application of an organic nitrogen source onto a soil was simulated. The output of that simulation, with known parameter values, was `sampled' to provide an error-free data set. The error-free data set was then manipulated to provide more realistic data by the addition of random error. These data were then used in the optimisations to estimate the mineralisation rate constant, k.
A successful optimisation may be defined as one in which the optimised parameter is both unbiased (where the optimised parameter has the same value as used in the original simulation) and precise (where the optimised parameter has an acceptably small confidence interval). If two or more parameters are simultaneously estimated then a further criterion is that the correlation between the optimised parameters should be small. Correlation between optimised parameters means that many combinations of the parameters could fit the data equally well, so casting doubt on the appropriateness of the optimised parameters.
Here we evaluate the precision and bias resulting from estimating k using measurements of nitrate concentrations alone, and the improvement resulting from adding a measurement of the cumulative nitrate flux or cumulative drainage obtained from a simple lysimeter. The usefulness of an additional sampling depth for nitrate concentration or a more intensive sampling regime is investigated. We also examine the problems associated with the need to optimise simultaneously a parameter controlling the water regime.
Modelling, scenario, and data sets
All computer modelling was done within APSIM, the Agricultural Productions systems Simulator (McCown et al. 1996). The modules of importance to the simulations described here are the water and solute movement module APSIM-SWIM (Huth et al. 1996) and a nitrogen transformation module, TWO-POOL. APSIM-SWIM, an adaptation of SWIMv2 (Ross et al. 1992; Verburg et al. 1996), employs a numerical implementation of Richard's equation for water movement and the advection-dispersion equation for solute movement. TWO-POOL is a nitrogen transformation module based on the model descriptions published by Johnsson et al. (1987) and Hutson and Wagenet (1992). The primary difference between TWO-POOL and the earlier models is that the faeces-N pool was omitted, and therefore there were only 2 pools of organic nitrogen, namely labile and resistant (corresponding to litter and humus in Johnsson et al. 1987).
The scenario used in the simulations was the application of sludge from an oxidation pond used to treat dairy waste to a bare soil. The sludge application provided 300 kg N/ha, with 90% of the nitrogen as labile organic N and the remainder as ammonium. The soil properties used were those determined for a freely draining Red Kandosol (Isbell 1996) at a site near Wagga Wagga, New South Wales (Bond et al. 1997). Meteorological data from the Wagga Wagga site were also used. Irrigation was supplied whenever a 16 mm soil water deficit accumulated, with sufficient water applied to provide a 20% leaching fraction. No plants were grown, and potential evaporation was taken as the measured pan evaporation with a pan coefficient (mm of potential evaporation per mm of evaporation from the pan) of 1.0 mm/mm. Actual evaporation was determined by the ability of the soil to meet that demand. Values of the parameters in the TWO-POOL model were as given by Johnsson et al. (1987) with the exception of the rate parameter for the mineralisation of the labile pool, k, which was set to 7.0x[10.sup.-2]/day. Short-term experiments are not suitable for determining the rate parameter for mineralisation from the resistant pool so that the parameter was set at a `default' value of 7.0x[10.sup.-5]/day, 3 orders of magnitude less than the parameter for the labile pool. After testing revealed almost no sensitivity of nitrate concentration at 1 m deep to the nitrification rate parameter, it was left at 0.2/day, the value recommended by Johnsson et al. (1987).
The above simulation was run for 182 days with the required outputs (namely, nitrate concentration in soil solution at depths of 0.5 and 1 m, cumulative nitrate flux past 1 m, and cumulative drainage past 1 m) reported at daily intervals. In practise, such data would be acquired by porous cup samplers or soil cores. The simulated nitrate concentrations are shown in Fig. 1 as the solid lines and are referred to as the underlying profile of nitrate concentration. Total drainage for the simulation period, D, and total nitrate leached, M, were also noted. `Observations' of nitrate concentration in the soil water at 1 m were obtained from the underlying profile at intervals of 2 weeks for the first 5 sampling times and then every 28 days until 182 days after the sludge application. These 10 samples are without error compared with the underlying profile, and are shown in Fig. 1 as closed circles. They will be referred to as the `no-error' data set. This rather sparse sampling strategy was chosen to represent the constraints to field sampling that normally apply when using distant experimental sites.
[Figure 1 ILLUSTRATION OMITTED]
For each sampling time, 5 additional concentrations containing `error' were calculated from the no-error data to represent replicates within the experimental site. Observations of nitrate concentration always include error arising from a variety of sources including spatial variability, the inability of some sampling devices to access all of the soil's water and solute, and analytical error. In the case of nitrate leaching, the former 2 sources usually dominate the total.
Many processes which cause variability in the field lead to skewed distributions (Parkin et al. 1988), but here we added normally distributed error to the `no-error' data set. Before data from skewed distributions are used for fitting parameters they must first be transformed so that the errors about the mean are normally distributed (Bates and Watts 1988). For that reason we chose to apply normally distributed error. Random, normally distributed noise was added to each of the 10 points in the `no-error' data set. The concentrations with error, C, were calculated from
(1) C = C'[1 + 0 [multiplied by] 5 N(0, 1)]
where C' is a concentration value from the `no-error' data set, and N(0, 1) represents a number drawn at random from a normally distributed population with a mean of 0 and a variance of 1. Random values that resulted in negative concentrations were discarded. The factor 0-5 controls the coefficient of variation of the resultant concentration values. The process of calculating C for each of the C' values (10 in total) was repeated 5 times to produce 5 replicate data sets with incorporated error. The coefficient of variation of the data was 46-5%, and was in the middle of the range for field data compiled by Parkin et al. (1988). The additional profiles, numbered 1-5, are shown in Fig. 1 as the open symbols of data set 2.
One of the scenarios investigated for improving the estimate of k was the measurement of nitrate concentrations at a second depth. To investigate this possibility, 5 replicate sets of nitrate concentrations at 0-5 m, with random noise, were obtained from the error-free simulation using the same technique as described above. These data are also shown in Fig. 1 as data set 3. Another scenario tested was more frequent sampling, for which a different set of 1 m concentrations was obtained from the error-free simulation results. In this case, the sampling frequency was weekly and resulted in 27 rather than 10 samples during the 182 days of the experiment. Five replicate data sets, totalling 135 observations, were created by the addition of random error as described above. These data are shown in Fig. 1 as data set 4.
In addition to the nitrate concentration profiles described above, 2 other measures were extracted from the simulation. These were M, the total mass of nitrate leached past 1 m, and D, the total drainage passing 1 m depth. The data values described above were assembled in various combinations into 6 data sets as summarised in Table 1.
Table 1. Summary of the six data sets used in the optimisations Data set Name Additional data 1 No error -- 2 With error -- 3 Extra sampling depth -- 4 Weekly sampling -- 5 Data set 2 plus measure M M 6 Data set 2 plus measure M and D M and D Data set No. of values 1 10 2 50 3 100 4 135 5 51 6 52
All parameter estimation was carried out with PEST (Watermark Computing 1994), a model-independent parameter estimation program that uses the Gauss-Marquardt-Levenberg method of non-linear optimisation. A uniform initial guess for k of 5.0x[10.sup.-2]/day and parameter bounds of 0.1 to 50x[10.sup.-2]/day were supplied for each optimisation.
Where replicate data are available, there are 2 methods with which to use them for optimisation. One method is to use all of the data at once to produce a single optimised value of k where no attempt is made to distinguish data from the different replicates. This will be referred to as the `massed data' method. For data set 2 (Table 1), it results in 50 concentration observations contributing to the calculation of the objective function. A second method is to use each replicate set of 10 concentration observations individually to produce 5 estimates of k, which are then combined to obtain a k appropriate for the entire data set. Since it was not obvious a priori which scheme is most appropriate, both were tested.
To test the second option, optimising to a single concentration profile at a time, the data must be assigned to 5 concentration profiles. This was done in 2 ways to represent the extremes of autocorrelation that might be found in data collected in the field. One extreme is where there is no serial correlation of the data. To achieve this, at each sampling time the concentrations were randomly assigned to 5 replicate profiles. This will be referred to as the `random data'. The other extreme is where a site that yields the highest concentration at 1 sampling time will always exhibit the highest nitrate concentration. This was mimicked by ranking the concentrations at each sampling time and using those ranks to assign data to the 5 replicate profiles. This will be referred to as the `ordered data'. The constant, k, was optimised using each of the replicates in the random and ordered data and the k for the data set, [k.sub.dataset], calculated by taking the median of the replicate k values, [k.sub.replicate]. The 95% confidence interval about [k.sub.dataset] was calculated from its variance, which was obtained by
(2) Var([k.sub.dataset]) = E[Var([k.sub.replicate])] + Var[E([k.sub.replicate])]
where Var is variance and E is the expected value (Dyson et al. 1990). The parameter values and their confidence intervals obtained for the random, ordered, and massed data were 9.6 [+ or -] 6 [multiplied by] 6, 8.14 [+ or -] 4.3, and 5.1 [+ or -] 1.9 x [10.sup.-2]/day. As there was no significant difference between the 3 methods tested, the massed technique was used in all further analysis as it was easiest to implement and avoided further difficulty of whether the data should be arranged in a random or ordered manner.
The first optimisations attempted to find the value of the mineralisation rate constant, k, from data sets 1 and 2 (Table 1). The purpose of these was to investigate the implications resulting from using concentration data subject to the degree of variability found naturally.
In order to address the difficulty of estimating k from data set 2, a revised experimental design might either include a second sampling depth or have an increased frequency of sampling. The first of these options was tested using data set 3 and the second by using data set 4.
The next set of optimisations investigated the utility of the measurement of the total amount of nitrate leached by adding M to the calculation of the objective function. The same optimisations as described above were repeated. The single value of M was given an arbitrary weight of 5, so that deviations of the simulated value from the `observed' value of M contributed 5 times as much to the total sums of squares as an individual concentration measurement. This was equivalent to regarding deviations in M (kg/ha) as being 5 times more important than a deviation of the same magnitude from a concentration observation (g/[m.sup.3]). This has implications for the required reliability of the techniques used to measure M, so the effect of the weight applied to M was investigated by repeating the optimisation using the combined data set with a weight of 1 or 2 applied to M to supplement the weights of 0 and 5 already assessed.
The fifth set of optimisations investigated the interaction between error in the measurement M and the weight applied to M. Here the error in M was assigned by arbitrarily subtracting 5% or 10% of the true value.
A final set of optimisations was carried out to investigate the effect of uncertainty in the soil water regime on parameter estimation. All of the previous optimisations assumed that the amount and temporal pattern of drainage was without error. Often it is necessary to optimise simultaneously 1 or more of the parameters controlling the simulation of the water balance in addition to k. The water balance affects nitrate leaching directly by changing the drainage, and therefore transport of nitrate to groundwater, and also indirectly by affecting the rate of mineralisation through soil water content and temperature. Evapotranspiration is often not directly measured in field experiments and it is necessary to obtain the pan coefficient as well as the mineralisation rate parameter from the data. To investigate the interaction between estimation of 2 parameters simultaneously, the pan coefficient, E, was optimised along with k using data set 6. The lower and upper bounds for E were 0.0 and 3.0 mm/mm and the initial value was 0.7 mm/mm. In the calculation of the objective function, there was no error in M and it was given a weight of 5. Combinations of weights of 0, 1, 2, or 5 on D with the effect of a 0, 5, or 10% underestimate of D on the optimisation were investigated.
For each optimisation, the root mean square error (RMSE) was computed using the deviations of predicted concentration from the appropriate no-error data set concentration. To allow comparison of optimisations performed using different data sets or data weighting, the RMSE was calculated using only the concentration data, even if M or D was used in the optimisation. Because the RMSE was calculated using concentration data from the no-error data set, it is not a measure of the goodness-of-fit of the optimisation, which had been given data with error. Rather it is a measure of how good the optimisation was at getting true concentration values despite the error in the data. This is necessary to have a basis for comparison of the results.
Results and discussion
The various data sets generated by the addition of random noise to the error-free simulation are shown in Fig. 1. Several statistics relating to the data sets containing error are given in Table 2, showing they had coefficients of variation slightly less than 50% and displayed a degree of variation consistent with that observed in field studies (Parkin et al. 1988; Lord and Shepherd 1993).
Table 2. Descriptive statistics for the data sets containing error M is the mass of nitrate beneath the concentration-cumulative drainage curve calculated by trapezoidal integration. The true values of M for the 0.5 and 1 m data are 297 and 288 kg/ha. CV, coefficient of variation
Data set CV of M from no- M from median of no. concn data error data data sets with error (%) (kg/ha) (kg/ha) 2 46 263 271 3 48 300 273 4 48 270 267
As a principal aim of modelling field studies such as those considered here is to simulate nitrate leaching, the amount of nitrate leached was calculated using data from the no error data and the median values of data sets 2, 3, and 4. With the relatively infrequent sampling, little improvement can be made on straight-line interpolation, or trapezoidal integration, between the data points. This interpolation can, however, lead to some error in the estimate of M, so here it is called apparent M, or M. The M calculated during simulation follows the underlying profile shown in Fig. 1, rather than being calculated from sparse points assuming a linear interpolation. Nevertheless it is useful to compare the trapezoidal integrations and this has been done in Table 2. It can be seen that M for the data sets with error is reasonably close to M for the error-free data, although there is a bias toward low M, particularly for the sparse 0.5 and 1 m data sets. If sampling ceases before the pulse of nitrate has completely passed a given sampling depth, the resulting low M may lead to an unrealistically low estimate of k. In our data sets, the nitrate concentrations had fallen to background levels for about 30 days before the last sample, so this bias is unlikely to directly cause low estimates of k.
Effect of error in the observations
The results of the first round of optimisations, where k was estimated from data sets with or without error, are given in Table 3. Also given in Table 3 are the 95% confidence intervals (CI), as a percentage of the estimated k, and the RMSE for each optimisation.
Table 3 Effect of three experimental options on the precision and bias of estimated k values
RMSE and the total number of measurements for the 5 replicate profiles (N) are also shown. The true value of k is 7.0 x [10.sup.-2]/day
Experimental k CI(A) RMSE N option ([10.sup.2]/day (%) Data set 2 5.1 37.3 4.1 50 Data set 3 6.0 22.5 53.5 100 Data set 4 6.9 18.8 18.8 135 Data set 5 6.9 2.9 0.1 50(B)
(A) 95% confidence interval expressed as a percentage of the parameter value.
(B) In addition, 1 measurement of the total amount of nitrate leached from the lysimeter is required
Where there was no error in the nitrate concentration observations, even with the infrequent sampling, the optimisation was able to produce an excellent estimate of k. This might be expected given that the integral beneath the underlying profile is not substantially different from the integral under the no-error observations. This is in agreement with the work of Lord and Shepherd (1993), who found that, for data with a single peak, as few as 4 samples were sufficient to calculate M within 10% of an intensively sampled data set.
When error was added to the data used for the optimisation, the estimated values of k showed considerable bias and a large confidence interval. A large value of CI should be regarded as an overestimate of the true value of CI because their calculation is based on linearity assumptions that might not apply to the full range of the interval (Watermark Computing 1994). In these cases, the CI should not be taken literally but should instead be regarded as indicating that the parameter is very poorly defined. The bias and lack of precision in k resulting from data set 2 indicates that, with data of this nature, 5 replicates are insufficient to define k adequately. In this example, random error was added to the samples. Many of the causes of variability and uncertainty in the field arise from sources where bias in the measurement may be introduced (Addiscott 1996). This bias is likely to be inherited by any parameter optimised using the data.
Evaluation of experimental strategies to improve the estimate of k
Strategies for providing data with better capability of producing unbiased and precise estimates of k are either to measure nitrate concentrations at a second depth (data set 3) or to increase the frequency of sampling (data set 4). The results of testing these options are given in Table 3. Either of the options improved both the precision and unbias of k. Sampling at an additional depth showed the least improvement in predicting k with respect to bias, whereas sampling at a weekly interval led to an estimate of k close to the true value. Both options resulted in CIs that were about 20% of the estimated value. While this is an improvement over data set 2, it is not sufficient to predict reliably mineralisation resulting from field applications of organic material.
Effect of including M in the objective function
The optimisations using data set 2 described above were repeated, this time including M in the calculation of the objective function (data set 5) with the application of an arbitrary weight of 5. The results of these optimisations are given in Table 3. Addition of M to the objective function significantly improved the estimate of k, which showed less bias and was more precisely defined than in the optimisations without M. This indicates the potential usefulness of simple lysimeter techniques, independent of the concentration data, that provide the mass of solute leached. However, in these optimisations, M was without error and weighted highly, indicating a high degree of confidence in the value. Before proceeding further, it is necessary to test both the effect of the weight on M and the effect of any error in the measurement of M.
Choice of the weight applied to M in the optimisations affects the estimate of k, as shown in Fig. 2. A weight of zero means that M does not affect the calculation of the objective function and, therefore, in that case data set 5 becomes effectively no different from data set 2. When a weight of 1 was applied, the estimate of k was closer to the true value, compared with a weight of zero, and the CI decreased from 37.3 to 3.8x[10.sup.-2]/day. Improvements in precision and bias of k were gained by increasing the weight on M further. For practical purposes, the unbias and precision in estimates of k resulting from a weight of 2 would be acceptable for simulation of processes in field experiments. This implies that any method used to obtain M in the field should be twice as reliable as a near-peak measurement of solute concentration. Here reliability is interpreted as the net result of analytical and technical capabilities and also spatial variability.
[Figure 2 ILLUSTRATION OMITTED]
The above parameter estimation assumed that the value of M had no error. In reality this is never the case. The next set of optimisations therefore examined the interaction between the error in the measurement of M and the weight applied to M in the optimisation. These results are summarised in Fig. 3. Not surprisingly, error in M led to significant error in the optimised value of k. However, the CIs on those estimates were similar at all error levels, and decreased markedly as the weighting on M increased. The optimisation for M containing 10% error and a weighting of 5 resulted in a value of k that was half the true value but was well defined with a small CI. This shows that the apparent precision, indicated by the CI, can be quite misleading. Thus, it is clear that if M is to be included in an optimisation then it must be an accurate measurement of the amount of solute leached. The technique used to measure M must be accurate, and as M is likely to be spatially variable, there must be sufficient replication to capture the true value of M. This may affect the usefulness of including M in the optimisations in many cases.
[Figure 3 ILLUSTRATION OMITTED]
Use of the value of M, such as might be measured using a lysimeter, resulted in a confidence interval of only 3%, compared with a 19% confidence interval for the weekly sampling (Table 3). Use of a lysimeter has the further advantage of being cheaper as it requires fewer samples to be analysed than additional nitrate concentration measurements and, important]y for distant research sites, significantly fewer sampling events. As discussed above, however, the accuracy and precision of the lysimeter option is strongly dependent on the lysimeter's ability to determine M accurately.
Effect o/uncertainty in simulation of the water regime
A set of optimisations were performed using data set 6 to investigate the interaction between estimation of k and simultaneous optimisation of the pan coefficient, E. The results of these optimisations are shown in Table 4. In the simulation used to generate the underlying profile, E was arbitrarily assigned a value of unity, and so a successful optimisation should return this value with a small CI. Additional considerations are that k should also be well defined and that any correlation between the 2 parameter estimates should be small.
Table 4. Estimates of k and E with their 95% confidence interval for data set 6
Also shown are the RMSE of the optimisation and the correlation between the optimised parameters. The true value of k is 7.0 x [10.sup.-2]/day, and of E is 1.0 mm/mm. See the text for details of the optimisations
Error Weight applied to D in D 0 1 2 5 k ([10.sup.-2]/day) 0% 6.9 6.9 6.9 6.9 -5% 10.1 10.1 10.1 -10% 35.6 39.2 39.5 E (mm/mm) 0% 1.0 1.0 1.0 1.0 -5% 1.2 1.2 1.2 -10% 1.5 1-5 1.5 RMSE 0% 0.1 0.1 0.1 0.1 -5% 2.8 2.9 2.9 -10% 4.3 4.3 4.4 Error Weight applied to D in D 0 1 2 5 % CI of k(A) 0% 17.4 4.3 3.6 2.9 -5% 3.0 2.5 2.0 -10% 1.0 0.6 0.5 % CI of E(%) 0% 28.6 5.8 2.9 1.2 -5% 4.8 2.4 1.0 -10% 4.1 2.1 0.8 Correlation of k and E (%) 0% 98.4 74.7 49.5 22.3 -5% -10%
(A) 95% confidence interval expressed as a percentage of the parameter value.
When D was not included in the calculation of the objective function, both k and E were identified with sufficient unbias but the precision of both parameters was unacceptable. The estimated value of E ranged between 0.7 and 1.3 mm/mm. This variation in E is sufficient to cause considerable differences in the simulated water regime. For example, if the soil were sufficiently wet that evaporation occurred at the potential rate, drainage at the lower end of the confidence interval would be about half that at the upper end of the confidence interval. In addition to the poor precision, the correlation between k and E was close to 100%, further indicating the lack of certainty in the identification of the parameters. Correlation between parameters found by simultaneous optimisation implies that several combinations of the 2 parameters would produce equally good fits to the optimisation data, thus leading to considerable uncertainty about the most appropriate parameter values. This correlation can cause difficulty when applying the parameters beyond the conditions to which they were originally defined.
When D without error was included in the optimisation, both k and E were successfully identified with acceptable unbias, even with a weight of only 1 applied. The principal change caused as the weight applied to D increased from 1 to 5 was that the correlation between k and E decreased from an unacceptably high value of 75% to an acceptably low value of 22%.
As the error in D increased the estimates of k and E became rather biased. For each level of error in D, changes in the weight applied to D primarily affected the confidence interval on the estimate, and had little effect on the parameter value. Interestingly the precision of the estimates increased as bias increased, meaning that caution should be applied to the maxim which states that a parameter with a small confidence interval is a sign of a good optimisation. From this it is clear that in order to find reliable estimates of mineralisation parameters it is essential that water movement is independently and accurately determined. Simultaneous monitoring of the movement of a tracer may aid in this.
If nitrogen transformation models are to be used with confidence for evaluating land management scenarios, then it is essential to collect site-specific information when determining the important parameters controlling mineralisation and leaching. It is possible, on theoretical and pragmatic grounds, to simplify most simulations to the point where only a few parameters need be optimised from experimental data, the rest being measured, independently estimated, or otherwise approximated. Experiments carried out with the aim of providing data to estimate rate parameters must be carefully designed if the parameters are to be identified with sufficient precision and unbias.
The analysis presented here shows that variability in nitrate concentrations results in considerable bias and lack of precision in k if infrequent measurements of nitrate concentration in the soil water alone are used. A single measurement of M during an experiment can improve the estimate of k, but only if it is accurately measured. In the scenario examined here, improvement in the value of k only occurred if the error in M was less than about 2%.
There are measurement techniques under development that would seem promising in this regard. The important characteristic of such a measurement is that it should sample the entire mass of solute leaving a defined volume of soil to allow accurate measurement of solute leached. Potentially useful methods include drainage lysimeters using intact soil, particularly with recent advances in their collection and installation (Moyer et al. 1996), and mini-lysimeters with ion exchange resin traps (e.g. Sakadevan et al. 1994). Care is needed in the interpretation of the data from mini-lysimeters because they are of shallow depth and often change the water regime. Work by Wyland and Jackson (1993) indicates that care in the design on the ion exchange trap is also essential.
The analysis presented here has also shown that the need to optimise a parameter describing the water regime simultaneously whilst estimating k adversely affects the unbias and precision of both estimates. Preference should be given to experimental designs that allow the parameters affecting the water balance simulation to be estimated independently of the optimisation of k. Here a properly designed lysimeter would again be a valuable addition to the experimental design. The results also suggest that optimisation of k by inverse methods will only be successful when other pathways of nitrogen loss, for example plant uptake and denitrification, are either quantified independently or known to be negligible.
Collection of nitrate concentration at a second depth, or more frequently from the one depth, improved the unbias and precision of the estimate of k, but neither strategy was as effective as an accurate measurement of M. The need for a short sampling interval has practical implications for experiments at distant research sites.
This work has also demonstrated the value of using model scenarios to aid the design of experiments and to estimate the bias and precision in parameter values obtained by optimisation when using field data.
The software for the modelling environment and some of the modules used in this study were developed by the Agricultural Production Systems Research Unit, a collaboration of CSIRO Tropical Agriculture and Queensland Department of Primary Industries. We wish to thank Mr John Hargreaves (CSIRO Tropical Agriculture) for assistance with programming the APSIM module interface, and Mr Neil Huth (CSIRO Tropical Agriculture) and Dr Kirsten Verburg (CSIRO Land and Water) for their work in incorporating SWIMv2 into APSIM.
This work was partially funded by the Land and Water Resources Research and Development Corporation under the project CWW17 `Australian Design Method for Sustainable Land Treatment of Rural Industry and Sewage Effluent'.
This work forms part of the CSIRO Wagga Wagga Effluent Plantation Project, which is partially funded by the Land and Water Resources Research and Development Corporation, the Murray Darling Basin Commission, and the NSW Department of Land and Water Conservation, and provided with in-kind support by Wagga Wagga City Council and Tahara Pastoral.
We thank Dr Dennis Rolston (University of California, Davis) and Dr Bronwyn Harch (CSIRO Mathematical and Information Sciences) for helpful comments on an earlier version of this manuscript.
Addiscott, T. M. (1996). That uncertain feeling: measuring and modelling nitrogen losses by leaching. Land and Water Resources Research Corporation, Occasional Paper No. 08/96, Canberra.
ARMCANZ and ANZECC (1995). Draft management guidelines for dairy sheds--December 1995. Commonwealth of Australia, National Water Quality Management Strategy.
Bates, D. M., and Watts, D. G. (1988). `Nonlinear Regression Analysis and Its Applications.' (Wiley: New York.)
Bond, W. J., Smith, C. J., and Ross, P. J. (1998). Field validation of a water and solute transport model for the unsaturated zone. In `Shallow Groundwater Systems: 2. Flow and Solute Transport Models'. (Heise Verlag: Hannover.)
Brandi-Dohrn, F. M., Dick, R. P., Hess, M., and Selker, J. S. (1996). Suction cup sampler bias in leaching characterization of an undisturbed field soil. Water Resources Research 32, 1173-82.
Cameron, K. C., Smith, N. P., McLay, C. D. A., Fraser, P. M., McPherson, R. J., Harrison, D. F., and Harbottle, P. (1992). Lysimeters without edge flow: An improved design and sampling procedure. Soil Science Society of America Journal 5 , 1625-8.
Chapman, S. L. (1996). Soil and solid poultry waste nutrient management and water quality. Poultry Science 75, 862-6.
Cook, M. G., Hunt, P. G., Stone, K. C., and Canterberry, J. H. (1996). Reducing diffuse pollution through implementation of agricultural best management practices--A case study. Water Science and Technology 33, 191-6.
DiStefano, F. J., and Gholz, H. L. (1986). A proposed use of ion exchange resins to measure nitrogen mineralization and nitrification in intact soil cores. Communications in Soil Science and Plant Analysis 17, 989-98.
Dyson, J. S., Jury, W. A., and Butters, G. L. (1990). The prediction and interpretation of chemical movement through porous media: The transfer function model approach. EN-6853 (Electric Power Research Institute: Palo Alto, CA.)
Helyar, K. R., and Porter, W. M. (1989). Soil acidification, its measurement and the processes involved. In `Soil Acidity and Plant Growth'. (Ed A. D. Robson.) pp. 61-101. (Academic Press: London.)
Huth, N. I., Keating, B. A., Bristow, K. L., Verburg, K., and Ross, P. J. (1996). SWIMv2 in APSIM: An integrated plant, soil water, and solute modelling framework. In `Proceedings 8th Australian Agronomy Conference, Toowoomba, 1996'. (Ed. M. Asghar.) p. 667. (Australian Society of Agronomy: Carlton, Vic.)
Hutson, J. L., and Wagenet, R. J. (1992). `LEACHM: Leaching Estimation and Chemistry Model: A Process Based Model of Water and Solute Movement Transformations, Plant Uptake, and Chemical Reactions in the Unsaturated Zone. Continuum Vol. 2. Version 2.' (Water Resources Institute, Cornell University: Ithaca, NY.)
Isbell, R. A. F. (1996). `The Australian Soil Classification.' Australian Soil and Land Survey Handbook Series No. 4. (CSIRO Publishing: Collingwood, Vic.)
Jarvis, S. C., Stockdale, E. A., Shepherd, M. A., and Powlson, D. S. (1996). Nitrogen mineralization in temperate agricultural soils: Processes and measurement. Advances in Agronomy 57, 187-235.
Johnsson, H., Bergstrom, L., Jansson, P-E., and Paustian, K. (1987). Simulated nitrogen dynamics and losses in a layered agricultural soil. Agriculture, Ecosystems and Environment 18, 333-56.
Larocque, M., and Banton, O. (1996). Using field data and simulation modelling to determine nitrogen cycle parameters. Soil Science Society of America Journal 60, 1840-5. Lord, E. I., and Shepherd, M. A. (1993). Developments in the use of porous ceramic cups for measuring nitrate leaching. Journal of Soil Science 44, 435-49.
McCown, R. L., Hammer, G. L., Hargreaves, J. N. G., Holzworth, D. P., and Freebairn, D. M. (1996). APSIM: A novel software system for model development, model testing, and simulation in agricultural systems research. Agricultural Systems 50, 255-71.
Metherell, A. K., Harding, L. A., Cole, C. V., and Parton, W. J. (1993). CENTURY soil organic matter model, technical documentation, agroecosystem Version 4.0. Great Plains Research Unit Technical Report No. 4. USDAARS, Fort Collins, CO.
Moyer, W. M., Saporito, L. S., and Janke, R. R. (1996) Design, construction, and installation of an intact soil core lysimeter. Agronomy Journal 88, 253-6.
Parkin, T. B., Meisinger, J. J., Chester, S. T., Starr, J. L., and Robinson, J. A. (1988). Evaluation of statistical estimation methods for lognormally distributed variables. Soil Science Society of America Journal 52, 323-9.
Popp, T., and Wichmann, W. (1996). Reducing nitrate in water--The European experience. Fertiliser and Lime Research Centre Occasional Report No. 9. pp. 31-43. Massey University, Palmerston North, New Zealand.
Raison, R. J., Connell, M. J., and Khanna, P. K. (1987). Methodology for studying fluxes of soil mineral-N in situ. Soil Biology and Biochemistry 19, 521-30.
Ross, P. J., Bristow, K. L., Bailey, S. W., and Smettem, K. R.J. (1992). Using SWIM (Version 2) to model soil water and solute movement. In `Agronomy Abstracts, Minneapolis, Minnesota, 1992'. (American Society of Agronomy: Madison, WI.)
Sakadevan, K., Hedley, M. J., and Mackay, A. D. (1994). An in situ mini lysimeter with a removable ion exchange resin trap for measuring nutrient losses by leaching from grazed pastures. Australian Journal of Soil Research 32, 1389-400.
Schnabel, R. R., Messier, S. R., and Purnell, R. F. (1993). An evaluation of anion exchange resin used to measure nitrate movement through soil. Communications in Soil Science and Plant Analysis 24, 863-79.
Verburg, K., Ross P. J., and Bristow K. L. (1996). SWIMv2.1 User Manual. CSIRO Aust. Division of Soils, Divisional Report No. 130.
Watermark Computing (1994). `PEST: Model-Independent Parameter Estimation.' (Watermark Computing: Brisbane.)
de Willigen, P., and Ehlert, P. A. I. (1996). Nitrogen and phosphorus balances and losses from the soil in the Netherlands. Fertiliser and Lime Research Centre Occasional Report No. 9, pp. 7-27. Massey University, Palmerston North, New Zealand.
Wyland, L. J., and Jackson, L. E. (1993). Evaluating nitrate recovery by ion-exchange resin bags. Soil Science Society of America Journal 57, 1208-11.
Yamaguchi, T., Moldrup, P., Rolston, D. E., and Hansen, J. Aa. (1992). A simple inverse model for estimating nitrogen reaction rates from soil column leaching experiments at steady water flow. Soil Science 154, 490-6.
Manuscript received 25 June 1997, accepted 19 September 1997
V. O. Snow(AB) and W. J. Bond(C)
(A) CSIRO Land and Water, PB 2, Glen Osmond, SA 5064, Australia.
(B) Centre for Groundwater Studies, PB 2, Glen Osmond, SA 5064, Australia.
(C) CSIRO Land and Water, PO Box 1666, Canberra, ACT 2601, Australia.
|Printer friendly Cite/link Email Feedback|
|Author:||Snow, V.O.; Bond, W.J.|
|Publication:||Australian Journal of Soil Research|
|Date:||Jan 1, 1998|
|Previous Article:||Recharge estimation for the Liverpool Plains.|
|Next Article:||Nitrogen fluxes in surface soils of 1-2-year-old eucalypt plantations in Tasmania.|