Assessing nitrogen fluxes from dairy farms using a modelling approach: A case study in the Moe River catchment, Victoria, Australia

Agricultural Water Management 178 (2016) 37–51

Contents lists available at ScienceDirect

Agricultural Water Management journal homepage: www.elsevier.com/locate/agwat

Assessing nitrogen fluxes from dairy farms using a modelling approach: A case study in the Moe River catchment, Victoria, Australia T. Thayalakumaran ∗ , A. Roberts 1 , C. Beverly, O. Vigiak 2 , N. Sorn, K. Stott Department of Economic Development, Job, Transport and Resources, 32 Lincoln Square Nth, Carlton, VIC 305, Australia

a r t i c l e

i n f o

Article history: Received 15 September 2015 Received in revised form 1 September 2016 Accepted 11 September 2016 Keywords: N loss pathways Waterway water quality Catchment heterogeneity Parameter sensitivity

a b s t r a c t Assessment of nitrogen (N) loss forms and pathways from farming systems is important for improved understanding of potential off-farm impacts on high value environmental assets. The objective of this study was to estimate N losses in different pathways in dairy systems across the range of climate, soil and farm management by using west Gippsland (Victoria, Australia) as case study area, and to characterise the sensitivity of the adopted model parameters. We combined the point scale models DairyMod and Howleaky to estimate dissolved N (DN) and particulate N (PN) loads in runoff, and N leaching (LN) in deep drainage from representative dairy farms in west Gippsland. Monte Carlo error propagation with Latin hypercube sampling was performed to identify sensitive model parameters and assess potential uncertainty in N load predictions. The combined model was capable of simulating climate-soil-animalpasture management interactions and estimating DN, PN and LN at an annual scale; which were estimated at up to 18 kg-N ha−1 , 15 kg-N ha−1 and 312 kg-N ha−1 , respectively. The combined model demonstrated that more intensive feeding as mixed ration, and nutrient budgeting that takes into account the fertiliser equivalent of recycled nutrients can achieve an increase in milk production by up to 13% and a decrease in N loads by up to 31% compared to the intensive system in the case study catchment. Soil type and farm management explained much of the variability (up to 76%) observed in LN and DN loads, whereas climate and soil type had significant influence on PN loads (62–77%). Year-to-year variation, particularly under dry conditions had a marked influence on N loads. Soil N, vegetation cover, rooting depth and soil maximum drainage rate must be well characterised in order to reduce potentially high uncertainty in the estimation of N losses in heterogeneous catchments. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Over the last two to three decades, as is the case in many countries, Australian dairying has intensified by increased stocking rates, greater reliance on supplementary feeding and increased nitrogen fertiliser usage (Gourley et al., 2012; Christie et al., 2012). This intensification not only increased milk production but also has resulted in greater on-farm nutrient surpluses contributing to increased degradation of water quality in waterways and water bodies and emissions of greenhouse gases (de Klein and Eckard, 2008; Bryant et al., 2011; Doole and Pannell, 2012; Smith and Western, 2013).

∗ Corresponding author. E-mail address: [email protected] (T. Thayalakumaran). 1 Present address: Natural Decisions Pty. Ltd., Melbourne, Australia. 2 Present address: European Commission, Joint Research Centre (JRC), Institute for Environment and Sustainability (IES), Via E. Fermi 2749, Ispra (VA), I-21027 Ispra, VA, Italy. http://dx.doi.org/10.1016/j.agwat.2016.09.008 0378-3774/© 2016 Elsevier B.V. All rights reserved.

While soil nutrient surplus is the “source factor” indicating the potential for nutrient losses, environmental conditions such as climate and landscape characteristics may facilitate or hinder the off-site transport of the surplus nutrients via various pathways (runoff, leaching and gaseous emissions). Estimating farm scale nutrient losses in different pathways and forms, and understanding the contribution of environmental conditions and farm management practices to loss pathways is essential both at small spatial scales (farms) and larger regional scales. At the farm scale, balanced on-farm management practices suited to the environmental conditions may be adopted to reduce losses and maximise production. At the regional or catchment level, this information can be used to inform and evaluate policy options, and help target locations for intervention or intensification to off-set nutrient contributions (Barns et al., 2009). Links between dairy intensification and N concentrations in waterways have been established (Smith et al., 2013) however, there exists limited knowledge about the environmental nitrogen impacts of dairy farms of varying complexity in Australia. A recent review concluded that insufficient research has been

38

T. Thayalakumaran et al. / Agricultural Water Management 178 (2016) 37–51

Nomenclature Climate zones FS Foot slopes HR High rainfall zone Soils I II III IV V VI VII

Friable earth Red iron-rich structured earth Heavy structured earth Texture contrast & gradational soils Heavy earthier and texture contrast soils Sodic texture contrast soils Cracking clays & wet soils

Periods Dry Wet

1990 to 1996 1997-2005

Dairy farm management systems System 1 Sys1 Sys2 System 2 Sys3 System 3 System 4 Sys4 Sys5 System 5 Model outputs DNL Dissolved nitrogen load LNL Leached nitrogen load Particulate nitrogen load PNL Q-runoff DD Deep drainage E Soil erosion Model parameters Air dry water content in soil AD CNb Curve number bare Curve number reduction at 80% cover CNr cv Response curvature Dissolved nitrogen concentration in water DNC DNf Multiplication factor for nitrate content in the surface soil DNSoil Nitrate content in the surface soil (0–2 cm) Water content in soil at field capacity FC GC Green cover GCf Green cover multiplication factor k Parameter for soil nitrogen and water nitrogen mixing Leaching efficiency LE LNC Leached nitrogen concentration in water LNf Multiplication factor for nitrate content in the soil layer contributing leaching TNSoil Total nitrogen content in the surface soil (0–2 cm) LNSoil Nitrate content in the soil layer contributing leaching Maximum drainage rate MD NERSoil Soil nitrogen enrichment ratio Residual cover RC RCf Residual cover multiplication factor RD Root depth RD Root depth multiplication factor Water content in soil at saturation Sat SW Soil water between air dry and saturated water content in soil

TNf WP

Multiplication factor for total nitrogen content in the surface soil Water content in soil at wilting point

conducted in the area of N leaching and some specific areas of surface N runoff in Australia, and recommended more research to quantify N leaching under different soil and climatic conditions representing Australian dairy production systems (Burkitt, 2014). Given the difficulties in constructing experiments to cover sufficient environmental heterogeneity, simulation models that integrate complex biological and physical processes are used to predict production benefits and nutrient losses from diverse agricultural systems, to assess best management practice impacts and to evaluate economic-environmental trade-offs (Cameira et al., 2007; Bryant et al., 2011; Vogeler et al., 2012). Nutrient export may vary greatly for combinations of climate, land use, soil type and farm management at the paddock/grid level. Despite the wide variety of available models (RZWQM (Ahuja et al., 2000), SGS (Johnson et al., 2003), APSIM (Keating et al., 2003), EPIC(Gassman et al., 2005), Howleaky (McClymont et al., 2007), DairyMod/EcoMod (Johnson et al., 2008), GLEAMS (Devereux, 2012), and IFSM (Rotz et al., 2012)), no models sufficiently simulate all the aspects of the dairy farming systems and their contribution to nutrient export in all pathways and forms. For example DairyMod and EcoMod can simulate paddock/farm scale complex pasture-animal-nutrient dynamics, predict marketable outputs, and N losses in leaching and gaseous forms (Eckard et al., 2006; Cullen et al., 2008; Snow et al., 2009), but cannot simulate N losses in runoff. In contrast, the Howleaky model (McClymont et al., 2007) simulates the effect of soil and land management on water balance, soil erosion and nutrient exports, but cannot simulate soil-pasture-animal-nutrient interactions. Understanding parameter sensitivity and its contribution to uncertainty in model estimates is key to making sound, robust environmental decisions and maximising the benefits attained from such decisions (Jakeman and Letcher, 2003; Payraudeau et al., 2007; Refsgaard et al., 2007). Monte Carlo simulations are commonly used to analyse how the distributions of input data and parameter errors are propagated through a simulation model, and how they affect model outputs (Balakrishanan et al., 2005; Miller et al., 2006). To minimise the number of Monte Carlo simulations, a Latin hypercube sampling (LHS) approach can be used which adopts an integrated random and stratified sampling strategy that considers all sections of the parameter space efficiently; thus it is efficient for dealing with many input variables and heterogeneous environmental characteristics (Vachaud and Chen, 2002; van der Keur et al., 2008). The objective of the study was to estimate N losses in different pathways in dairy systems across the range of climate, soil and farm management using a modelling framework with application to a west Gippsland (Victoria, Australia) catchment as a case study. DairyMod was combined with the Howleaky model to estimate N losses in the dominant pathways and forms. Model parameter sensitivity associated with water and N pathways and fluxes was assessed with a global sensitivity analysis procedure. 2. Materials and methods 2.1. Study area The Moe River catchment, with an area of 577 km2 within the Gippsland Lakes basin, south east Australia (Fig. 1) was selected for this study. The Moe River, a tributary of the Latrobe River has been

T. Thayalakumaran et al. / Agricultural Water Management 178 (2016) 37–51

39

Fig. 1. Major soil types (a) and land use categories (b) in the Moe River catchment.

identified as a major contributor of nutrient loads to the Gippsland Lakes (Hancock et al., 2007) which is an iconic environmental asset and Ramsar Site. Total nitrogen (TN) contribution from the Moe River to the Gippsland Lakes has been estimated as 124 t yr−1 for the period between 1975 and 2005 (Vigiak et al., 2013). Two agro-climatic zones exist in the catchment: the foot slopes (FS, at altitude of 500–1000 m a.s.l.; average rainfall of 750 mm yr−1 ) and a high rainfall zone (HR, at altitude of 1000–1500 m a.s.l.; average rainfall of 1050 mm yr−1 ). Precipitation occurs in all months, with the least in February and most in September and October (www.bom.gov.au). The upland of the catchment is dominated by friable earths; the western slopes by iron rich structured earth (25% of catchment); the north-eastern slopes by texture contrast and some gradational soils (43%); whereas the southern slopes are dominated by texture contrast, sodic soils (10%). The main plain is dominated by vertosols and hydrosols (19%) (Vigiak et al., 2013; Fig. 1b). Dairy (30%) and beef/dairy mix (44%) are the main enterprises of the catchment, and the rest of the catchment is occupied by forest plantations, horticulture and urban settlements (Fig. 1b). Land management of dry land dairy farms was classified into four main management systems representing the range in fertiliser usage, supplementary feed usage, stocking rates, and milk production observed in the catchment (Table 1). Details of method development for defining dairy management systems and the characteristics are in Stott et al. (2013).

2.2. Models This study adopted arguably the most commonly used models in Australia to simulate dairy enterprises (DairyMod) and runoff processes (Howleaky), importantly both models adopt daily time-steps and a cascading bucket approach to represent soil water redistribution. DairyMod (Johnson et al., 2008) is a one-dimensional (vertical), biophysical model developed for New Zealand and Australian grazing farming systems. The model simulates pasture growth and utilisation by grazing animals, milk production, water and nutrient dynamics in soils and N leaching for a range of pasture, animal, fertiliser, and irrigation management options. A number of paddocks are integrated into a farm by controlling the grazing of animals around the paddocks, with supplements made and/or bought in, depending on pasture quantity and quality relative to animal feed demand. The model is described in detail elsewhere (Johnson et al., 2008). Howleaky (McClymont et al., 2007) is a one-dimensional (vertical), field-scale model designed to simulate the impact of different land uses and management on water balance, soil erosion and water quality. The model uses the water balance and erosion components of PERFECT v3 (Littleboy et al., 1999). A recently included N module has had limited testing on pastoral systems (Thayalakumaran et al., 2013), but the conceptualisation for estimating field-scale N losses in runoff and leaching are similar to other contemporary models such as SWAT, EPIC and GLEAMS.

40

T. Thayalakumaran et al. / Agricultural Water Management 178 (2016) 37–51

Table 1 Attributes of four dry land dairy management systems of increasing management intensity in Moe River catchment. Note that the characteristics are defined for the milking area. System

Intensive (Sys1) Moderately intensive (Sys2) Moderately extensive (Sys3) Extensive (Sys4)

N fertilisation

Supplement fed

Stocking rate

Milk production

kg ha−1 yr−1

t DM cow−1 yr−1

cows ha−1

L cow−1 yr−1

L ha−1 yr−1

2.2 1.9 1.6 1.5

3.0 2.5 2.2 1.9

6300 6000 5200 4200

18,900 15,200 11,500 7,900

a

300 200b 100c 50d

a = April, May, June, July, August, September, October and November; b = April, May, August, September, October and November; c = May, September and October; d = May and September.

The modelling approach adopted in this study were as follows 1. Adapt a calibrated DairyMod developed for the Gippsland region by Eckard et al. (2006) 2. Conduct a series of simulations representing variations in climate, soil and farm management using DairyMod 3. Collate the time series estimates of pasture growth, topsoil nutrient concentrations and subsoil nutrient concentrations from each of the DairyMod simulations. 4. Incorporate the DairyMod outputs collated time series estimates from each of the DairyMod simulations into Howleaky to predict time series of nitrogen loads in various pathways. A schematic of the approach linking DairyMod to Howleaky in order to represent all nutrient loss pathways is shown in Fig. 2. The specific model component interactions are detailed as follows. For easy readability, the model parameters and intermediate outputs are in italics and the main model outputs are in regular text. In Howleaky, dissolved N concentrations (DNC , mg l−1 ) and loads (DNL , kg-N ha−1 ) in runoff are calculated as: DNC (t) = DNsoil (t) × k × [1 − exp(−c v × Q (t))]

(1)

DNL (t) = DNC (t) × Q (t)

(2)

(mg-N kg−1 ) is the nitrate concentration in the surface

where DNsoil soil layer (0–2 cm) contributing to runoff; k is a parameter used to describe mixing of soil N with runoff water; cv describes the response curvature and Q is runoff (mm). Particulate N loads (PNL , kg-N ha−1 ) in runoff is calculated as: PNL (t) = E(t) × TNsoil (t) × NERsoil × 10−6

(3)

(kg ha−1 );

where; E = erosion TNsoil = total nitrogen content of topsoil (mg-N kg−1 ); NERs oil = soil Nitrogen Enrichment Ratio. Nitrate concentration (LNC , mg l−1 ) and loads (LNL , in kg-N ha−1 ) in leaching is estimated as a function of nitrate concentration in water in the lower soil layer, and deep drainage (DD in mm). Nitrate concentration in soil water is calculated as, LNC (t) =

BNsoil (t) SW (t)

LNL (t) = LNC (t) × DD(t) × LE

(4) (5)

where BNsoil (mg-N kg−1 ) is daily (or monthly) nitrate concentration in the soil layer contributing to leaching (if known) or the 0–10 cm soil depth; SW (mm) is the volume of soil water between air dry water content and saturated water content of the entire soil profile or of the lower soil layer; LE is a leaching efficiency parameter, which partitions the soil water nitrate concentration into runoff and deep drainage pathways (ranges between 0 and 1); and DD is daily (or monthly) deep drainage (mm). In this application DairyMod generated the vegetation growth time series and soil N concentrations (DNsoil , TNsoil and BNsoil ) at various depths on a daily time step for all climate-soil-farm management combinations. DairyMod daily time series outputs were

then imported into Howleaky as monthly average time series. It means that each day of the month had the same soil N concentration. Howleaky simulations were then undertaken to estimate N losses in runoff and leaching as well as to assess model parameter sensitivity. In acknowledgement that soil N concentration and vegetation cover are influenced by water balance, the predictions of evapotranspiration, drainage and runoff estimates derived from the two models were verified for consistency. 2.3. Simulations In total 56 combinations of climate zone (2) × soil (7) × farm management (4) were modelled for 16 years (1990–2005) following a 10 year initialisation period (1980–1989). The 16 year simulation period was split into two hydrologically distinctive periods based on rainfall-runoff time series analysis: a relatively wet period (1990–1996), and a following dry (1997–2005). 2.3.1. Climate and soil data DairyMod requires daily rainfall, minimum and maximum temperature, minimum and maximum relative humidity, net radiation and vapour pressure. Howleaky climate data requirements are daily rainfall, pan evaporation, temperature and solar radiation. All the climate data for were extracted from the SILO website as daily time series (SILO, Australian Bureau of Meteorology- http://reg.bom.gov. au/). In the simulation period, FS received a mean annual rainfall of 838 mm yr−1 in wet period and 684 in the dry period (751 for 1990–2005), whereas in the HR climatic zone mean annual precipitation was 1182 mm yr−1 in wet period and 949 mm in dry period (1051 mm yr−1 for 1990–2005). The soil parameters were set according to the national database of soil properties (McKenzie et al., 2000; Table 2). DairyMod allows the soil profile to be discretised into 4 layers while three layers are recommended in Howleaky. The upper two soil layers were equal in thickness in the two models, but the deeper layer 3 in Howleaky was partitioned into layer 3 and 4 in DairyMod. Soil characteristics were the same in both models except for cracking clays and wet soils. Cracking clays and wet soils exhibit crack flow which was represented by assigning a maximum crack infiltration of 100 mm in Howleaky, while DairyMod cannot model this process. 2.3.2. Implementation of dairy farm management in DairyMod An existing calibrated DairyMod simulation developed for the Gippsland area using experiment data by Eckard et al. (2006) was adapted in this study and with modification to the attributes/parameters relevant to the Moe River catchment dairy systems. In all dairy management systems, the farm size was set to 110 ha with 28 paddocks. Rotational grazing of pasture (perennial rye grass and white clover mix) based on leaf stage was assumed for all paddocks. Cereal grains (barley, 12.3 MJ kg-DM−1 ) and pasture silage (9.4 MJ kg-DM−1 ) were used as concentrate and forage respectively for all dairy management systems. N content was set to 2.06% for barley and 2.6% for pasture silage. Dung and urine

T. Thayalakumaran et al. / Agricultural Water Management 178 (2016) 37–51

41

Fig. 2. A schematic of the approach linking DairyMod to Howleaky in order to represent all nutrient loss pathways at the farm scale.

removed from the milking area per day was set to 15% to reflect the time cows spent in laneways and in dairy sheds/yards (Aarons and Gourley, 2012). The equivalent N amount was removed from the milking area. Model parameters pertaining to the nitrogen cycle were as in Eckard et al. (2006). Measured organic carbon content in the dairy farms soils of Moe River catchment by Vigiak et al. (2012) were set as initial values. P fertilisation was set to 20 kg ha.yr−1 . To implement the 4 dairy management systems and to derive the system characteristics as in Table 1, mature animal body weight and peak milk production (L day−1 ) were set. Minimum and maximum daily concentrate and forage (kg animal−1 ) were varied to achieve the total amount of supplements fed (t DM cow−1 yr−1 ) of the four dairy management systems, while maintaining a concentrate: forage ratio of 3:1 for N. Nitrogen fertiliser (urea) was applied in several splits of 40 kg ha−1 or less depending on the dairy management systems. There were eight, five, three, and two applications per year for systems 1, 2, 3 and 4 respectively. Applications were on the 1st of the month in each year. 2.3.3. Generation of pasture cover time series for Howleaky In Howleaky, daily plant growth is simulated using either a static or dynamic approach. Because the effect of animal grazing on pasture growth could not be considered in Howleaky, the static approach was used for this study. The static approach required temporal patterns of green cover (GC, 0–100%), residue cover (RC, 0–100%), and root depth (RD). DairyMod simulated the daily time series of live shoot, dead shoot, litter and root biomass. Live shoot and root biomass were converted into daily time series of GC and RD respectively, and dead shoot combined with litter was converted to RC. A nonlinear relationship of cover (%) defined by 100 × 1.0 e0.6 × Leaf Area Index was used to convert dry mass into percent cover. Daily time series were then summarised into monthly average time series. The resulting monthly average time series were not significantly different between climate zones, soils or dairy management systems. Thus the same monthly average time series of GC, RC and RD were used in Howleaky for all systems. 2.3.4. Generation of soil N time series for Howleaky To simulate DNL , LNL and PNL Howleaky required time-series of DNsoil , BNsoil and TNsoil (Eqs. (1) to (5)). These were derived from DairyMod simulations: daily soil nitrate content in 0–2 cm was used as DNsoil ; soil nitrate concentration in the last 10 cm soil depths was used as BNsoil ; and total N content in the soil surface (0–2 cm) was used for TNsoil . Model interface limitations constrain

DairyMod to report only one user-defined paddock soil N content outputs. However, preliminary analysis showed that the variability in soil N content was not significantly different between the 28 simulated paddocks of the representative farm. Therefore, the DairyMod paddock for generating the daily time-series was randomly chosen within the simulated farm configuration. DairyMod daily time series of DNSoil , BNsoil and TNsoil were summarised into monthly averages. Monthly average DNsoil , BNsoil and TNsoil time series differed among the 56 climate-soil-system combinations. Thus, separate monthly average time series of DNsoil , BNsoil and TNsoil were imported in Howleaky for each combination. The TNsoil were checked with the measured TN values for consistency. Measured TNsoil varied from 6100 to 8300 mg kg−1 in Soil II & Soil VII, while all the other soils measured 4500 to 5800 mg kg−1 (Vigiak et al., 2012). It is acknowledged that the DNsoil in 0–2 cm soil might be variable however, there was no measured data available for validation. 2.3.5. Statistical analysis To estimate the effects of climate, farm management and soils, and their interaction on the model outputs an Analysis of Variance (ANOVA) was performed on all model outputs using Genstat14th edition. Variance for each of the stratum (climate, farm management, soils, climate.soils, climate.farm management, soils.farm management, climate.farm management.soils) and the total variance were used to estimate the percentage of variance explained. 2.4. Sensitivity analysis The Howleaky configuration adopted in this study has more than 20 model input parameters that describe soil water, vegetation and nutrient processes. In order to limit the computational effort to undertake a model sensitivity analysis a Latin Hypercubic Sampling (LHS) scheme was implemented in the SENSAN package of PEST (Doherty et al., 1995). This procedure allowed the assessment of the sensitivity of water and N flux estimates to Howleaky model parameters in all 56 climate-soil-system combinations. Following main steps were involved in the sensitivity analysis. 2.4.1. Defining model parameter boundaries and distribution functions A preliminary sensitivity analysis, using one-parameter-at-atime (OAT) approach guided the selection of sensitive model

42

T. Thayalakumaran et al. / Agricultural Water Management 178 (2016) 37–51

Table 2 Model parameters for the soils in the Moe River catchment. Model parameters in bold are common for both DairyMod and Howleaky models, non-bolded parameters are only applicable to the Howleaky model. Soil Groups

Friable earth

Red iron-rich Structured earth

Heavy structured Texture contrast & Heavy earthier andSodic texture earth gradational soils texture contrast contrast soils soils

Cracking clays & wet soils

Soil codes a Field texture Three layer depths (mm)

I L 100 300 1350 6 8 10 23 23 25 45 45 41 63 63 53 100 100 100 0.9–1.2 11 80 25 19.4 8.5 3.75 0.4 1

II vFSCL 100 300 500 7 9.5 12 23 24.5 26 45 43 41 64 58.5 53 75 75 50 0.9–1.2 2.4 85 25 11.5 9.1 4 0.25 1

III CL 100 300 600 3 4.5 6 23 26 26 48 41 41 63 51 51 50 50 25 1.2–1.5 4.3 85 30 30.7 9.4 4 0.3 1

VII LC 100 300 1200 7 8.5 10 31 31 31 45 44 43 52 49 46 5 5 1 1.2–1.4 2.8 75 20 1 9 4 0.18 1

Airdry water content (AD)

Wilting point (WP)

Field capacity (FC)

Saturated water content (Sat)

Maximum drainage rate (MD) (mm day−1 ) Bulk density b Organic carbon (%) CN-bare CN-reduction 80% cover Slope (%) Stage1 SoilEvap U (mm) Stage2 SoilEvapCona (mm) USLE K USLE P a b

IV ZCL 100 300 950 3 7.5 12 9 19.5 30 22 27.5 33 33 34 35 25 25 5 1.6–1.7 0.5 90 30 20 6.75 3.5 0.55 1

V FSCL 100 300 750 3 3 12 10 10 7 27.5 27.5 21 35 35 35 25 25 5 1.6–1.7 3.9 90 30 2 8 3.6 0.4 1

VI ZCL 100 300 1200 2 6.5 11 9 19.5 30 22 27.5 33 34 34 34 5 5 1 1.6–1.7 4.0 90 30 1.1 8.6 3.8 0.25 1

L = loam; vFSCL = very fine sand clay loam; CL = clay loam; ZCL = silt clay loam; FSCL = fine sand clay loam. Vigiak et al. (2012).

parameters. In total, 19 model parameters were identified as sensitive for water and N fluxes. The parameters were grouped as: - Soil N parameters: DNsoil , BNsoil , TNsoil , NERsoil , cv and k. - Vegetation parameters: GC, RC and RD - Soil hydraulic properties for the three soil layers: Sat, FC, WP, MD, CNb and CNr (Table 2). All parameters were treated as independent to each other, as no observation data was available to estimate correlations between parameters. The initial probability distribution functions (pdf) and statistical boundaries of the 19 parameters were based on the soil database, DairyMod outputs, and expert judgement. For soil parameters (Table 3), the 5th and the 95th percentile values in the state wide soil coverage data base by McKenzie et al. (2000) and local soil pit data were used. Howleaky required some rules in soil water retention parameters (Sat, FC, WP) to be respected: FC cannot exceed Sat, and WP cannot exceed FC. To preserve these rules, the distribution functions (pdfs) were generated for Sat, for the ratio of Sat: FC, and for the ratio of FC:WP for the soil layers. From the distributions of these three meta-parameters, FC and WP pdfs were derived. Vegetation cover and soil N parameters in Howleaky consisted of monthly average time series. Because it was too cumbersome to generate distributions for each month of the time series, multiplication factors were used to reflect the variability in applying a shift (upwards or downwards) of the entire time-series. GCf, RCf, RDf, DNf, BNf and TNf were the factors used for GC, RC, RD, DNsoil , BNsoil and TNsoil time series respectively. For all factors a uniform distribution was assumed, but the boundaries of each factor were different and were set based on the DairyMod daily outputs for the 56 combinations (Table 4). Factors varied for the wet and the

dry periods. The distributions of soil N parameters in Table 4 were not varied with climate, soil, or management systems, as it was assumed that variations in soil N parameters due to these factors were captured in the base time series of DNsoil , BNsoil and TNsoil for the 56 combinations. van der Keur et al. (2008) recommended that the number of LHS parameter sets was either four-thirds the number of parameters or two to five times the number of parameters involved. We generated 50 parameter sets for the 19 parameters which satisfied the first recommendation and 2 times the number of parameter number. Note for all the 56 combinations of simulations 50 parameter sets were generated. 2.4.2. Quantification of sensitivity in model outputs The 50 parameter sets generated by LHS were propagated to obtain model outputs via Howleaky-SENSAN interface. Therefore, for each of the 50 parameter sets, 16 years of model outputs (Q, DD, DNL , PNL and LNL ), were generated. The model outputs from the reference base parameter set (single value) were compared with the model outputs from the 50 parameter sets (distribution) to assess their agreement. That is, we would expect that the model outputs from the reference base parameter set to be close to the median. Furthermore, non-parametric MannWhitney test was performed using Genstat 14th edition to verify whether two independent distributions have the same centre, that is the median. The Standardised Ranked Regression Coefficients (SRRC) method was used to further assess model sensitivity (van der Keur et al., 2008). With this method, a regression technique is used to rank the model parameters’ sensitivity based on the relative magnitude of the regression coefficients, which measure the influence of each parameter on the model output. Due to differences in units and magnitudes of the parameters and of model outputs, all model out-

T. Thayalakumaran et al. / Agricultural Water Management 178 (2016) 37–51

43

Table 3 Soil parameter statistics and probability density functions (pdf) for the sensitivity analysis.

pdf Soil statistics mean I SD lower upper

MD1 log

MD2 log

MD3 log

Sat1 normal

Sat2 normal

Sat3 normal

Sat1-FC1 uniform

Sat2-FC2 uniform

Sat3-FC3 uniform

FC1-WP1 uniform

FC2-WP2 uniform

FC3-WP3 uniform

CNb normal

CNr normal

300 3.16 – –

300 3.16 – –

300 3.16 – –

63 8 – –

63 8 – –

53 6.5 – –

0.18 – 0.13 0.23

0.18 – 0.13 0.23

0.12 – 0.1 0.16

0.22 – 0.16 0.28

0.22 – 0.16 0.28

0.16 – 0.13 0.21

80 1.7 – –

25 1.7 – –

II

mean SD lower upper

300 3.16 – –

300 3.16 – –

300 3.16 – –

64 8 – –

58.5 6.5 – –

53 6.5 – –

0.19 – 0.14 0.24

0.16 – 0.12 0.19

0.12 – 0.1 0.16

0.22 – 0.16 0.27

0.19 – 0.14 0.23

0.15 – 0.13 0.20

80 1.7 – –

25 1.7 – –

III

mean SD lower upper

100 3.16 – –

100 3.16 – –

30 3.16 – –

63 7.8 – –

51 4.5 – –

51 4.5 – –

0.15 – 0.11 0.18

0.1 – 0.08 0.11

0.1 – 0.08 0.11

0.25

0.15

0.15

0.18 0.3

0.11 0.17

0.11 0.17

80 1.7 – –

30 1.7 – –

IV

mean SD lower upper

30 3.16 – –

30 3.16 – –

3 3.16 – –

33 5.5 – –

34 5.5 – –

35 5 – –

0.11 – 0.07 0.14

0.07 – 0.04 0.08

0.02 – 0.01 0.03

0.13 – 0.08 0.17

0.08 – 0.05 0.1

0.03 – 0.02 0.04

90 1.7 – –

30 1.7 – –

V

mean SD lower upper

30 3.16 – –

30 3.16 – –

3 3.16 – –

– –

– –

– –

0.13 – 0.07 0.16

0.08 – 0.04 0.09

0.02 – 0.02 0.05

0.13 – 0.07 0.17

0.08 – 0.07 0.17

0.03 – 0.03 0.14

90 1.7 – –

30 1.7 – –

VI

mean SD lower upper

10 3.16 – –

10 3.16 – –

3 3.16 – –

34 3.75 – –

34 3.75 – –

34 4.75 – –

0.12 – 0.09 0.16

0.07 – 0.05 0.08

0.01 – 0.008 0.012

0.13 – 0.1 0.18

0.08 – 0.06 0.1

0.03 – 0.02 0.04

90 1.7 – –

30 1.7 – –

VII

mean SD lower upper

5 3.16 – –

5 3.16 – –

1 3.16 – –

52 5.3 – –

49 5 – –

46 5 – –

0.07 – 0.06 0.09

0.05 – 0.04 0.06

0.03 – 0.02 0.04

0.14 – 0.13 0.17

0.13 – 0.11 0.16

0.12 – 0.1 0.15

75 1.7 – –

20 1.7 – –

MD: maximum drainage rate (mm day−1 ); Sat: saturated water content (mm); Sat-FC: ratio of saturated water content to field capacity (unit less); FC-WP: ratio of field capacity to wilting point (unit less); CNb: CN number for bare soil (%); CNr: CN number reduction for 80% cover (%). Numbers next to the text indicate soil layers. SD means standard deviation. Table 4 Nitrogen (N) and vegetation parameter boundaries for the sensitivity analysis. A uniform probability density function was assumed for all the parameters, and the parameters are unit less. default

lower

upper

N parameters DNf BNf TNf cv k NERsoil

1 1 1 0.2 0.5 1

0.1 0.3 0.5 0.1 0.1 1

2.7 2.3 1.9 0.8 1.0 2

Vegetation parameters GCf RCf RDf

1 1 1

0.7 0.7 0.7

1.2 1.2 1.2

TNf: factor for total nitrogen content of topsoil; BNf: factor for nitrate concentration in the soil layer contributing to leaching; DNf: factor for the nitrate concentration in the surface soil layer; NERsoil: soil nitrogen enrichment ratio;k:parameter used to describe mixing of soil N with runoff water; cv: describes the response curvature of runoff losses; GCf: factor for green cover; RCf: factor for residue cover; RDf: factor for root depth.

puts and parameters were first normalised. Then a multiple linear regression was applied on the rank-transformed data. The derived regression coefficients are reported as the SRRC indices. 3. Results 3.1. Water and soil losses A comparison of predicted runoff and deep drainage for various climate and soil combinations is provided in Table 5. Modelled mean annual runoff ranged from zero to 314 mm, being greatest on poorly draining soils (VI and VII), and least in well-draining soils (I, II and III). Runoff was greater in the HR zone compared

to the FS, particularly in poorly draining soils. A clear distinction in runoff generation between the dry and wet periods was also observed, with the wet period having three times the runoff as a percentage of rainfall compared with the dry period. The impact of dry-wet period on runoff generation was greater in HR than in FS, particularly in poorly draining soils. Soil erosion ranged from 0 to 3.8 t ha−1 yr−1 ; and broadly followed the runoff variation, being low in well-draining soils and high in poorly draining soils. Annual runoff to annual erosion in the dairy paddocks of Moe River catchment can be generalised by, Erosion (t ha−1 ) = 0.0107 RO (mm) + 0.0106 with an R 2 of 0.86

44

T. Thayalakumaran et al. / Agricultural Water Management 178 (2016) 37–51

Fig. 3. Box plots of annual N leaching on foot slopes, FS and high rainfall zones, HR (D = dry and W = wet) for various soils in Sys1 and Sys2. Note the scales on the Y axis.

Table 5 Mean annual runoff and deep drainage for the period of 1990–2005 (wet = 1990–1996; dry = 1997–2005). Numbers in brackets are the standard errors of means. Soils

Runoff (mm)

FS HR

Deep drainage (mm)

FS HR

I

II

III

IV

V

VI

VII

wet dry wet dry

1 (0.4) 0 1 (0.3) 0.2 (0.1)

1(0.4) 0 1(0.3) 0.2 (0.1)

4(1) 0 2(0.6) 1 (0.4)

31(6) 11(2) 49(7) 24(4)

30(6) 9(2) 31(5) 15(3)

93(9) 43(6) 314(16) 130(11)

91(7) 38(5) 297(15) 126(10)

wet dry wet dry

117(11) 26(6) 435(13) 222(13)

127(12) 33(7) 444(13) 231(13)

126(11) 34(7) 445(13) 231(13)

155(9) 73(8) 460(13) 261(12)

99(9) 30(6) 420(13) 218(12)

76(5) 30(4) 173(5) 139(5)

50(6) 8(2) 152(6) 112(5)

FS: foot slopes; HR: high rainfall zones.

Mean annual deep drainage ranged from 0 to 460 mm, with greater drainage in HR than in FS (Table 5). As expected, the deep drainage was the greatest in well to moderately draining soils, and least in soil VII. Unlike the runoff, the deep drainage was more impacted by dry and wet periods in FS than in the HR, particularly in well to moderately draining soils (IV and V).

3.2. Losses in runoff and leaching Mean annual dissolved N loads up to 17.9 kg-N ha−1 yr−1 were estimated (Table 6). DN load estimates for I, II, III and Sys4 are not presented as they were negligible (<0.2 kg-N ha−1 ). DN loads were more impacted by the variation between the wet and dry periods (by as much as 3 times) than by climatic zone (Table 6). Regardless of the climatic zone or farm management systems, the DN loads were greater in poor draining soils due to high runoff. Intensively managed systems on poorly draining soils generated the maximum DN loads. The mean annual DN load was less than 1 kg-N ha−1 for any farm management systems on well-draining soils and for the extensive systems on moderately drainable soils.

Particulate N loads were low (<4 kg-N ha−1 yr−1 ) from all soils except for the moderately draining soils that produced up to 4.2 and 15.1 kg-N ha−1 yr−1 due to large amount of erosion however, the PN loads were not different between systems (data not shown). Mean annual leached N ranged from 0 to 312 kg-N ha−1 with a wide year-to-year variation (Fig. 3). N leaching increased with farm management intensification; mean annual LN averaged over soils and climate zones was the greatest in Sys1 (160 kg-N ha−1 yr−1 ) and nearly halved in Sys2 (70 kg-N ha−1 yr−1 ), further reduced by 5 fold in Sys3, and was the least in Sys4 (7 kg-N ha−1 yr−1 ). In intensive systems Sys1 and Sys2, LN was greater in HR than in FS (except in soil III and IV in wet periods). The absence of such a trend in less intensive systems (data not presented) indicates that soil N content rather than water was the limiting factor of N leaching. Soils with a low field capacity and moderate maximum drainage rates (III, IV, and VI) leached more N than the soils with high field capacity and maximum drainage rates (I and II). The soil VII leached the least amount of N due to very low maximum drainage and high field capacity (Fig. 3 and Table 2).

T. Thayalakumaran et al. / Agricultural Water Management 178 (2016) 37–51

3%

Sys1

Sys2

45

4%

43%

38% 58%

1%

2% 0%

0%

3%

Sys3

51%

Sys4

3%

18% 0%

LN

3%

DN

12%

0% 5%

PN Vol Den

76%

80%

Fig. 4. Average annual N loads portioned into leaching (LN), dissolved N in runoff (DN), particulate N in runoff (PN), volatilisation (Vol) and denitrification (Den) from Sys1 Sys2, Sys3 and Sys4.

Table 6 Estimated mean annual dissolved N loads in runoff in kg-N ha−1 for the period of 1990–2005 (wet = 1990–1996; dry = 1997–2005) from three dairy management systems and four soils. Climate

FS

Systems

Sys1

HR FS

Sys2

HR FS

Sys3

HR

Period

dry wet dry wet dry wet dry wet dry wet dry wet

Soil Moderately draining soils

Poorly draining soils

IV

V

VI

VII

2.7 6.7 3.5 4.7 1.9 4.6 2.4 3.2 0.9 1.6 1.2 1.4

2.4 6.6 2.7 2.8 1.8 4.0 1.7 2.0 0.8 1.7 0.7 0.8

7.8 12.1 9.9 14.4 6.1 8.8 6.4 9.3 4.1 3.5 3.1 2.8

9.6 17.0 12.6 17.9 5.4 7.5 7.5 11.0 2.7 1.8 3.3 2.1

FS: foot slopes; HR: high rainfall zones.

3.3. Gaseous N losses from DairyMod The DairyMod model estimated up to 28 kg-N ha−1 yr−1 as direct denitrification fluxes with the greatest rate for management system Sys1 on poorly draining soils, where high soil nitrate-N and water logging conditions facilitated the denitrification process. Volatilisation flux varied between management systems: it was greatest in Sys1 (62 kg-N ha−1 yr−1 ) and least in Sys4 (27 kg-N ha−1 yr−1 ). No differences were predicted between climatic zones or soils due to model setting. Based on the Eckard et al. (2006) DairyMod, 15% of the applied urea fertiliser and 30% of urine were lost to volatilisation at 20 ◦ C; these percentages are solely linked to the farm management systems and not to soil or other climate factors.

207 kg-N ha−1 yr−1 for the Sys1, Sys2, Sys3 and Sys4 respectively. On average 41%, 32%, 20% and18% of the total N inputs to the systems were lost in various pathways and forms from the Sys1, Sys2, Sys3 and Sys4 respectively. As shown in Fig. 4, nitrate leaching was the major N loss pathway to water ways from all farm management systems. The greatest percentage of total N input that was returned to the soil as urine and dung was observed in Sys 3 and 4 (both 61%) compared to the more intensive Sys 1 and Sys 2 (47% and 54%). These increased amounts of N returns to soils in the form of urine in Sys 3 & 4 led to greater percentage of total N losses via volatilisation (Fig. 4).

3.4. Total N inputs to exports

3.5. Effect of climate, soil and farm management on water and N fluxes

On average, the total N input (fertiliser N, supplementary feed, N fixed) to the farm management systems were 553, 400, 271 and

Climate, soil and farm management including the interactions, and the year to year variation explained up to 96% of the variation

46

T. Thayalakumaran et al. / Agricultural Water Management 178 (2016) 37–51

in water fluxes (Q, DD), and up to 87% of the variation in N fluxes (DN, PN and LN) (Fig. 5). The residuals reported in Fig. 5 indicate the unexplained variances. Results suggest that DD was the best explained and DN was the least explained model output. Combined influence of climate, soil and farm management including the interactions was greater in the wet period than in the dry period for all model outputs except for DN, for which the impact was slightly less in the wet period. Impact of interactions alone on the model outputs ranged from 4 to 19%. Year-to-year variation had greater impact than the climate zone on most model outputs (up to 16%) except for DD for which the reverse was true, i.e. DD was more affected by climate zone than the year-to-year rainfall pattern. The contribution of year-to-year variation was much greater in the dry period than in the wet period. Climate and soil contributed markedly to the variation of water fluxes with 69–88% of the explained variation in runoff, 75–89% in deep drainage. The effects were significant (P < 0.001). Soil was the dominant driver of runoff (up to 60% of the explained variation), while climate contributed to much of the variation in deep drainage (up to 62%). Soil and management system were the main factors explaining LN and DN loads. Conversely, the variation in PN loads was mostly explained by soil and climate (62–77%) all at P < 0.001 significance (Fig. 5b). 3.6. Sensitivity analysis results 3.6.1. Comparison of base simulation with stochastic simulation model outputs Figs. 6 and 7 present the LHS based Monte Carlo simulation means and coefficient of variation (CV%) in the model outputs. As shown by the CV%, the variations associated with N fluxes were greater than in their relative water fluxes, i.e. variation in LN was greater than in DD, and variations in DN and PN loads were greater than in Q. The variation in all model outputs was always greater in the dry period than in the wet period. Comparison between the means of the base simulation (Table 5) and the Monte Carlo simulation showed relatively small difference for deep drainage in well-draining soils (<20% difference) in both dry and wet periods. Conversely poor draining soils (V and VI) exhibited large difference of DD, up to 50% (Fig. 7). Small differences between means for well-draining soils and large difference for poor draining soils were also observed for LN (Figs. 5 a and 7). Difference in Q means was relatively large compared to the DD in all soils (up to 60%), and significantly greater in the dry period (Table 5 and Fig. 6).

Fig. 5. Percentage contribution to explained variation in annual (a) water fluxes as runoff (Q) and deep drainage (DD), and (b) N fluxes in leaching (LN), dissolved N in runoff (DN) and particulate N in runoff (PN) due to climate, soil, system, year and their interactions.

gen concentration in deep layer) was a very sensitive (ranked 1–4) and significant parameter for LN in all soils except VII. Soil hydraulic parameters, particularly the MD at the 1st and 3rd layers were significantly sensitive in moderate to poorly draining soils (IV, V, VI and VII). DN loads were very sensitive to the N parameters k and DNf (the soil N concentration in top soil). PN loads were sensitive to a mix of soil parameters and to the TNf in moderate to poorly draining soils (data not shown). 4. Discussion

3.6.2. Model parameter ranking Tables 7–9 present the estimated ranks for the parameters based on the absolute SRRC values. Only parameters that had significant SRRC values (P < 0.001) at least in one soil type are presented for clarity. For most soils except soil IV, GC (green vegetation cover) was an important and significantly sensitive parameter for runoff estimation (Table 7). Top soil hydraulic parameters (MD1, WP1) ranked high for soils with moderate to poor drainage capacity (IV, V, VI and VII), but the relationship was not always significant. Insignificant SRRC values means that the evidence for a relationship between the independent (MD1, WP1) and dependent variable is weak. All vegetation parameters (GCf, RCf and RDf) were very sensitive for DD in many soil types in both wet and dry periods (Table 8). Medium to high ranking (1–6) indicate that the RDf and GCf were among the most sensitive parameters for deep drainage in all soils except for soil VII. MD was the most important and significant parameter controlling deep drainage, particularly in moderate to poor draining soils (IV, V, VI and VII). As shown by the high ranks, LN was sensitive to vegetation parameters (GCf and RDf, Table 9) in most soils. Also BNf (soil nitro-

4.1. Predicted deep drainage The equivalent recharge rates based on the deep drainage estimates summarised in Table 5, when allowance is made for subsurface partitioning as described by Rassam and Littleboy (2003), are within the bounds defined by groundwater hydrograph responses and as reported by GHD (2010) and Beverly et al. (2014). 4.2. Predicted N fluxes The N fluxes predicted with the combined DairyMod-Howleaky model were generally consistent with N losses reported in literature in comparable dairy systems and environmental conditions. The annual total N export of 0–21 kg-N ha−1 estimated in this study is comparable to 0.06–12 kg-N ha−1 measured by Barlow et al. (2007) for the same soil (V) with similar farm management (Sys 2) characteristics in the Gippsland region, Victoria. DN was the major form of N loss in runoff with up to 84% of the total N as DN is also consistent with Barlow et al. (2007) who reported 88%.

T. Thayalakumaran et al. / Agricultural Water Management 178 (2016) 37–51

47

250 III

IV

V

VI

VII

Q (mm), DN and PN (kg-N ha-1)

15

200

150

18 23

100

50 50

0

35

31

78 43

89

72 10

75

dry

98

wet

98 78

68 82

85 84

dry

79

88

32

89 70 67

wet

Q

55

67 70 66

dry

wet

DN

PN

DD (mm) and LN (kg-N ha-1)

Fig. 6. Results of Monte Carlo simulated runoff (Q), dissolved N in runoff (DN) and particulate N in runoff (PN) as annual means (on the axis) and CV % (as numbers). Bars are not shown for soil I and II because of very low values.

500 74 I

400

II

III

IV

V

75

VI

300 36 52

200 61

100 98 51 80

34

47

23

79

59

84 69

77

55 58

71 91

34 27 31 98

0 wet

dry

dry

DD

wet LN

Fig. 7. Results of Monte Carlo simulated deep drainage (DD) and N leaching (LN) as annual means (on the axis) and CV % (as numbers). Bars are not shown for VII because of very low value.

Table 7 Runoff parameter ranking (based on the SRRC) in different soils and periods. Values in bold indicate the coefficients are significant (P < 0.001). Blank space means the ranking is higher than 10 (i.e. the parameter is not sensitive).

CNb CNr GCf MD1 MD3 RCf RDf WP1

Soil

I

Period

dry

wet

II dry

wet

III dry

wet

7

8 2 7

9 8 7

10 8 1 7

10

2 9

1 2 5 8

6 8 4

1 10 6

2 5 3

6

6

IV

1 9

7

dry

10 4 3 8 5 1

V wet

9 5 1 6 7 3

VI

VII

dry

wet

dry

wet

7

7

8

10

9 1

6 1 4

5 3 7 6 9 1

2 7 1 6 8 3

10 2

3

dry

wet

7 5

6 2

6

7 9 4

2

GCf: factor for green cover; RCf: factor for residue cover; RDf: factor for root depth; MD: Maximum drainage rate; WP: Wilting point; CNb: CN number for bare soil; CNr: CN number for reduction 80% cover. Numbers next to the text indicate soil layers. SRRC means standardised ranked regression coefficients.

Nitrogen leaching estimated in this study for farm management Sys2 on poorly draining soils in HR zone (93 kg-N ha−1 yr−1 ) was greater than the amount reported by Eckard et al. (2004) (up to 22 kg-N ha−1 yr−1 ) for a comparable dairy system on dermosolic, redoxic hydrosol in HR zone of the Moe River catchment. The reason for the very small amount of N leaching in Eckard et al. (2004) is unclear however. Conversely, LN estimated in this study

is comparable to those found in literature for similar dairy systems in other areas. The Ledgard et al. (1999) measured LN of up to 101 kg-N ha−1 yr−1 for a dairy system with a stocking rate of 3.3 cows ha−1 and 200 kg N fertiliser application on a free draining soil in New Zealand. Ledgard et al. (1999) system is comparable to farm management Sys1 (3.0 stocking rate, 300 kg-N ha−1 fertiliser) in HR zone, but with a more moderate fertiliser application

48

T. Thayalakumaran et al. / Agricultural Water Management 178 (2016) 37–51

Table 8 Parameter ranking (based on the SRRC) for deep drainage for different soil types in wet and dry periods. Values in bold indicate the coefficients are significant (P < 0.001). Blank space means the ranking is higher than 10. Soil

I

Period

dry

wet

10

6

CNb FC1 FC2 FC3 GCf MD1 MD3 RCf RDf Sat1 WP1 WP2

9 2 6 3 1 4 7

II

2 5 3 1 4 10

dry 3 7 10 4

8 6 9 1 2

III

IV

V

wet

dry

wet

dry

wet

10 3 8

10 6 5

7 4

10 4 2

9 4 2

4

2

2

6 10 8

6

5 9 3 1

6 10 3 1

7 6 9 2 1

3 1 4 7 8

9 3 1 6 8

7

VI

dry

wet

7 8 1 4 6 5 9 3 2 8

7 12 2 4 6 5 10 3 1 8

VII

dry

wet

dry

wet

7 10 1

7 8 3

5 7

5 15 6

5 4 3

4 6 1

3 1

6 2

5 2 9

3 1 10 8 8 2

8

8 4 2 10

GCf: factor for green cover; RCf: factor for residue cover; RDf: factor for root depth; MD: maximum drainage rate; Sat: saturated water content; FC: field capacity; WP: wilting point; CNb: CN number for bare soil. Numbers next to the text indicate soil layers. SRRC means standardised ranked regression coefficients.

Table 9 Parameter ranking (based on the SRRC) for N leaching for different soil types in wet and dry periods. Values in bold indicate the coefficients are significant (P < 0.001). Blank space means the ranking is higher than 10.

BNf CNb FC2 FC3 GCf MD1 MD3 RCf RDf Sat1 WP1 WP2

Soil

I

Period

dry

wet

dry

wet

dry

wet

dry

wet

3 7

1 5

2

1

2

1

3

3

10

2

9 10 6

8

9 3

5

3

2 6 10 9

2 5 10 10

8 7 5 1

9 7 6 1

4 1 10

II

4 2

10 7 8 4 1

III

10 7 9 6 2

6 1 4 9 8

IV

3

4 2 5 9 6

V dry

9 1 4 6 8 3 2 5

VI wet

9 1 4 6 7 10 3 2 8

VII

dry

wet

dry

wet

4 6 2

3 6 9

8

10

7 3 5

8 5 2

9 1

10 1

5 9

4

3 1

3 1

7 2 9

4 2 7

BNf: factor for nitrate concentration in the soil layer contributing to leaching; GCf: factor for green cover; RCf: factor for residue cover; RDf: factor for root depth; MD: maximum drainage rate; Sat: saturated water content; FC: field capacity; WP: wilting point; CNb: CN number for bare soil. Numbers next to the text indicate soil layer. SRRC means standardised ranked regression coefficients.

(200 kg-N ha−1 ), which resulted in LN of 128 kg-N ha−1 yr−1 in our study. LN of 7 kg-N ha−1 yr−1 estimated for the extensive system Sys4 in our study is below the 13 kg-N ha−1 reported by Betteridge et al. (2007) for dairy with 60 kg-N ha−1 application in Taupo region, New Zealand. Less LN in our study could be related to reduced fertiliser usage (50 kg-N ha−1 yr−1 ) and drainage amount compared to Betteridge et al. (2007). Temporal pattern in LN are consistent with the literature, in particular to Eckard et al. (2006) findings. The average annual denitrification loss of 10–19 kg-N ha−1 and volatilisation loss of 48 kg-N ha−1 yr−1 estimated in this study for farm management Sys2 in poorly draining soils are close to 13 kg-N ha−1 yr−1 and 57 kg-N ha−1 yr−1 respectively, reported by Eckard et al. (2003) for similar dairy system on hydrosol (poorly draining soil) in HR zone of the Moe River catchment. Other literature from New Zealand and Australia on N losses (Ledgard, 2014; Burkitt, 2014; Drewry et al., 2006) could not be compared to our study due to different environmental/system conditions or to insufficient information regarding these factors. Nevertheless, these literature comparisons though limited to a few data sets, indicate that the N fluxes predicted in this study are sensible. The combined model seems thus capable of capturing the important climate-soil-animal-pasture management interactions and to generate sensible N fluxes at least on an annual scale.

4.3. Effect of climate, soil and farm management on N fluxes Climate, soil and farm management significantly contributed to N fluxes, but to a different extent depending on N pathways and forms. Soil was the main driver of runoff and soil erosion, whereas climate explained much of the variation in deep drainage. It is noteworthy that one monthly time series of each GC, RC and RD was used in Howleaky for all the simulations regardless of climate, soil and farm management differences. This is because DairyMOD outputs of pasture ground cover were not sensitive to climate and soil type. Thus, the contribution of climate and soil type on N fluxes was primarily due to the absolute rainfall amounts and soil parameterisation. Soil type and farm management were the main contributors to LN and DN fluxes, but the variation in PN loads was mostly explained by soils and climate. Farm management was represented by monthly time series of soil N content (DNsoil , TNsoil and BNsoil ) in Howleaky. The time series of TNsoil were not significantly different between climate, soil or farm management types, thus PN losses were mainly driven by the runoff and soil erosion, hence by soil and climate. On the other hand, BNsoil and DNsoil were very different between farm management and soil types, hence much of the variation in LN and DN fluxes were explained by both farm management and soils. Greater contribution of the farm management

T. Thayalakumaran et al. / Agricultural Water Management 178 (2016) 37–51

than soils to LN and vice versa to DN may be related to the different time scales of the DN (daily) and LN (monthly or greater) generation processes compared to the variation of DNsoil and BNsoil . Daily or sub-daily time scale of the DN generation process also explains why DN loads were the model output least explained by the environmental conditions and farm management (i.e. 60%). The fact that climate was a dominant contributor to DD but not so much to LN suggests that the coincidence of DD with the availability of BNsoil is more important for controlling LN than the annual amount of DD. The effect of climate, soil, farm management, and year-to year variation on runoff, deep drainage, DN, PN and LN were better explained in the wet period than in the dry period. In the dry period we observed accumulation of soil N in DairyMod simulations, which was flushed during storm events. Varying accumulation of soil N and soil water between years may have caused the less explainable nature of the climate, soil and farm management on the outputs (Q, DD, DN, PN and LN) in dry period. This confuses the contributions of the primary drivers (climate, soil and farm management) on the water and N fluxes. The significant contributions of climate, soil and farm management to N fluxes highlight the importance of appropriately choosing management practices that are appropriate to environmental conditions and that a spectrum of management choices can be used to target locations for minimizing environmental impacts.

4.4. Sensitivity of N fluxes to model parameters Among the 19 model parameters included in the analysis, BNf was the most influential parameter of LN followed by the RDf and GCf in most soils. Among the soil hydraulic properties, MD was a sensitive parameter for DD and LN in soils IV, V, VI and VII, but less so for the I, II and III. These results were expected, because soils I, II and III had negligible (<4 mm) runoff and the soils have very limited influence on the water balance. Therefore, a single set of central parameter values (mean or median) can be used for these soils. However, in less water conductive soils (IV, V, VI and VII), the lower MD values become the limiting factor of infiltration process. DNf and k in most soils, and WP1 in moderate/poor draining soils (IV, V, VI and VII) were the sensitive parameters for estimating DN loads whereas the PN loads were very sensitive to TNf and GCf. The importance of vegetation parameters on N losses found in this error propagation analysis implies that the N fluxes are very sensitive to the climate. All these results are consistent with the variance analysis in Section 3.5 which found soil and climate as the main drivers of N transport variables (Q and DD), soil and farm management as the factors driving DN losses, and climate and farm management influencing PN loads. The significant error contribution of the vegetation, nitrogen and soil hydraulic parameters on the simulated DN, PN and LN, as shown by their large CV% indicates potential for considerable uncertainty in model outputs. The parameters that must be best characterised are the most sensitive ones. The sensitivity analysis results imply that there is less need to account for variability in soil hydraulic

49

parameters within I, II and III soils. Hence a small number of field measurements may be adequate to represent these soil types when estimating the mean annual N fluxes. Conversely, for IV, V, VI and VII great care is needed to account for the full variability in soil hydraulic parameters to obtain reliable predictions of DD and LN on annual basis. Variability in vegetation parameters becomes crucial only in dry years in the Moe River catchment. This is because the pasture ground cover in dairy farms has mostly been above 90% in the catchment. The variability of DNsoil , BNsoil and TNsoil needs to be fully accounted for if models of this type are to be used for estimating N fluxes. In these cases, measurements conducted in the study catchment rather than estimated from the national data base may improve the reliability of N output predictions.

4.5. Dairy system with reduced N losses and increased milk production Given that estimates of DN, PN and LN loads at an annual scale for a range of climate, soil and management combinations were sensible, the combined model was used to explore dairy systems that could produce the same or greater amount of milk production with less N loads compared to the intensive system Sys1. A hypothetical farming system (Sys5) was built based on the less common farms with greatest milk production, moderate fertiliser usage and high supplements usage observed in the Moe River catchment, and considering the farming practices and attributes that could be implemented in the combined model (Stott et al., 2013). Sys5 was constructed by modifying Sys1 by i) using bought-in mixed rations with high Metabolisable Energy (MJ kg-DM−1 ), 16–18% crude protein and 30–35% Neutral Detergent Fibre, and ii) using nutrient budgets that take account of the fertiliser equivalent of nutrients recycled in dung and urine (through use of a fully-concreted feed pad and effluent system upgrade). The total N applied was 190 kg ha−1 yr−1 ; 140 kg as N fertiliser plus 50 kg N equivalent of effluent. These were accommodated in DairyMod by increasing the dung and urine removed to 25% with the introduction of a feed pad. It was assumed that 50% of N volatilises while in effluent storage ponds. Effluent application was on the total milking area. Boughtin concentrates were increased to 2.8 t DM cow−1 yr−1 from 2.2 t DM cow−1 yr−1 in Sys1. Other characteristics remained the same as in Sys1. Table 10 shows the improvements in Sys5 compared to Sys1. In Sys5 the total N loads were reduced by 31% and 11% in FS and HR compared to Sys1, with a gain in milk production by 5 and 13% in FS and HR respectively. Increased predicted milk production (L ha−1 ) and decreased N loads ha−1 resulted in increased N efficiency compared to Sys 1 by 37% and 23% in FS and HR, respectively. The largest impacts were in N loads via leaching, with 43 to 68% reductions in FS and HR respectively. These results suggest that Sys5 is a potential dairy system for the Moe River catchment to achieve production benefits and improved environmental benefits compared with Sys1. Higher milk production and savings in chemical fertiliser may increase the net private benefits to the farmer,

Table 10 Comparison between dairy management Sys 1 and Sys 5 for annual milk production and N losses. Sys1

Milk production (L ha−1 ) Fertiliser application (kg-N ha−1 ) Animal N intake (kg-N ha−1 ) Total N losses (kg-N ha−1 ) N leaching (kg-N ha−1 ) N efficiency (g N loss L milk−1 ) FS: foot slopes and HR: high rainfall zones.

Sys5

FS

HR

FS

HR

18,963 300 465 242 160 12.7

19,044 300 489 273 194 14.3

20,055 190 521 166 90 8

21,561 190 564 242 62 11

50

T. Thayalakumaran et al. / Agricultural Water Management 178 (2016) 37–51

but a thorough cost benefit analysis that takes into account the increased investment in construction of a feed pad and effluent pond infrastructure is required. Notwithstanding these improvements, Vigiak et al. (2013) showed that even with 100% of the Sys5 in the Moe River catchment, the N loads would still be 40% greater than the current composition with 10% of Sys 1 (most intensive), 30% of Sys 2, 40% of Sys 3, and 20% of Sys 4 (least intensive). Overall, environmental impacts from intensification appear to only be partially mitigated by current technologies. 5. Conclusions This study has demonstrated the utility of combining DairyMod with Howleaky for exploring N fluxes in the Moe River catchment dairy farming systems. Despite the limited amount of published data available for model validation, the combined model was found capable of capturing the climate, soil, animal, pasture management interactions, and of estimating dissolved, particulate and leaching N loads at an annual scale. Further, the relative contribution/influence of situational and management factors in Moe River catchment to various N loss pathways and forms was illustrated by the model. From the modelling perspective, global sensitivity analysis undertaken in the study identified soil nitrogen concentration, vegetation cover (particularly the root depth) and maximum drainage rate as the most sensitive parameters for estimating annual N losses for the combined model. Thus, the analysis indicated where and under what climatic conditions these sensitive parameters need to be sufficiently parameterised to reduce errors in N loads estimates. Acquisition of locally measured data and sensitive model parameters that are specific to the study catchment may further enhance the validation of the models’ ability to simulate processes and reduce uncertainty in the estimated N loads. Using the combined model, though difficult to implement all possible management practices, it is predicted that increased use of purchased feed concentrates as a mixed ration, nutrient budgeting that takes into account the fertiliser equivalent of recycled nutrients, and use of a feed pad and additional effluent storage capacity, can achieve increase in milk production by up to 13% and decrease in N loads by up to 31% compared to a current intensive system in this Moe river catchment. Changes in net private benefits of the alternative dairy systems or management practices including infrastructure investments were not included and would need to be evaluated as part of decision-making. Overall, environmental impacts from intensification appear to only be partially mitigated by current technologies. Acknowledgements This work was funded by the Victorian Department of Primary Industries. We thank David Rees (Victorian Department of Primary Industries) for his support in characterising soils. The authors also thank the reviewers and the editors, whose comments and advice helped to improve the manuscript substantially. References Aarons, S.R., Gourley, C.J.P., 2012. Sustainable management of nutrient returns in excreta on grazed dairy soils. In: Burkitt, L.L., Sparrow, L.A. (Eds.), Proceedings of the 5th Joint Australian and New Zealand Soil Science Conference: Soil solutions for diverse landscapes. Hobart, pp. 461. Ahuja, L.R., Rojas, K.W., Hanson, J.D., Shaffer, M.J., Ma, L. (Eds.), 2000. Water Resources Publ., LLC, Highlands Ranch, CO., p. 372. Balakrishanan, S., Roy, A., Ierapetritou, M.G., Flach, G.P., Georgopoulos, P.G., 2005. A comparative assessment of efficient uncertainty analysis techniques for environmental fate and transport models: application to the FACT model. J. Hydrol. 307, 204–218. Barlow, K., Nash, D., Hannah, M., 2007. The effect of fertiliser and grazing on nitrogen export in surface runoff from rainfed and irrigated pastures in south-eastern Australia. Nutr. Cycl. Agroecosyst. 77, 69–82.

Barns, A.P., Willock, J., Hall, C., Toma, L., 2009. Farmer perspectives and practices regarding water pollution control programmes in Scotland. Agric. Water Manag. 96, 1715–1722. Betteridge, K., Hoogendoorn, C.J., Thorrold, B.S., Costall, D.A., Ledgard, S.F., Park-Ng, Z.A., Theobald, P.W., 2007. Nitrate leaching and productivity of some farming options in the Lake Taupo catchment. N. Z. Grassl. Assoc. 69, 123–129. Beverly, C., Hocking, M., Cheng, X., O’Neil, C., Schroers, R., Baker, S., 2014. Conceptualisation and calibration of a Gippsland basin groundwater model, Department of Sustainability and Primary Industries, Melbourne. Tech. Rep., 176. Bryant, J.R., Snow, V.O., Cichota, R., Jolly, B.H., 2011. The effect of situational variability in climate and soil choice of animal type and N fertilisation level on nitrogen leaching from pastoral systems around Lake Taupo, New Zealand. Agric. Syst. 104, 271–280. Burkitt, L.L., 2014. A review of nitrogen losses due to leaching and surface runoff under intensive pasture management in Australia. Soil Res. 52, 621–636, http://dx.doi.org/10.1071/SR13351. Cameira, M.R., Fernando, R.M., Ahuja, L.R., Ma, L., 2007. Using RZWQM to simulate the fate of nitrogen in field soil-crop environment in the Mediterranean region. Agric. Water Manag. 90, 121–136. Christie, K.M., Gourley, J.P., Rawnsley, R.P., Eckard, R.J., Awty, I.A., 2012. Whole-farm systems analysis of Australian dairy farm greenhouse gas emissions. Anim. Prod. Sci. 52, 998–1011. Cullen, B.R., Eckard, R.J., Callow, M.N., Johnson, I.R., Chapman, D.F., Rawnsley, R.P., Garcia, S.C., White, T., Snow, V.O., 2008. Simulating pasture growth rates in Australian and New Zealand grazing systems. Aust. J. Agric. Res. 59, 761–768. de Klein, C.A.M., Eckard, R.J., 2008. Targeted technologies for nitrous oxide abatement from animal agriculture. Aust. J. Exp. Agric. 48, 14–20. Devereux, J., 2012. Managing Litter Re-use for Minimal Nutrient Run-off to Surface Water RIRDC Publication No. 11/160, ISBN 978-1-74254-338-3. ISSN 1440–6845. Doherty, J., Brebber, L., Whyte, P., 1995. PEST: Model independent parameter estimation. Australian Centre for Tropical Freshwater Research James Cook University, Townsville. Doole, G.J., Pannell, D.J., 2012. Empirical evaluation of nonpoint pollution policies under agent heterogeneity: regulating intensive dairy production in the Waikato region of New Zealand. Aust. J. Agric. Res. Econ. 56, 82–101. Drewry, J.J., Newham, L.T.H., Greene, R.S.B., Jakeman, A.J., Croke, B.F.W., 2006. A review of nitrogen and phosphorus export to waterways: context for catchment modelling. Mar. Freshw. Res. 57 (8), 757–774. Eckard, R.J., Chen, D., White, R.E., Chapman, D.F., 2003. Gaseous nitrogen loss from temperate perennial grass and clover dairy pastures in south-eastern Australia. Aust. J. Agric. Res. 54, 561–570. Eckard, R.J., White, R.E., Edis, R., Smith, A., Chapman, D.F., 2004. Nitrate leaching from temperate perennial pastures grazed by dairy cows in south-eastern Australia. Aust. J. Agric. Res. 55, 911–920. Eckard, R.J., Johnson, I., Chapman, D.F., 2006. Modelling nitrous oxide abatement strategies in intensive pasture systems. Int. Congr. Ser. 1293, 76–85. GHD, 2010. West Gippsland CMA Groundwater Model, Transient model development report. Produced by GHD (Gutteridge Haskins & Davey) on behalf of the Victorian Government Department of Sustainability and Environment, Melbourne, Australia. https://ensym.dse.vic.gov.au/docs/WestGippsland TransientModelReport FINAL.pdf/ (accessed 22.08.16.). Gassman, P.W., Williams, J.R., Benson, W., Izaurralde, R.C., Hauck, L.M., Jones, C.A., Atwood, J.D., Kiniry, J.R., Flowers, J.D., 2005. Historical Development and Applications of the EPIC and APEX Models. Working 49 Paper 05-WP 397, Centre for Agricultural and Rural Development, Iowa State University, Ames, Iowa 50011-1070. www.card.iastate.edu. Gourley, J.P., Dougherty, W.J., Weaver, D.M., Aarons, S.R., Awty, I.A., Gibson, D.M., Hannah, M.C., Smith, A.P., Peverill, K.I., 2012. Farm-scale nitrogen, phosphorus, potassium and sulphur balances and use efficiencies on Australian dairy farms. Anim. Prod. Sci. 52, 929–944. Hancock, G., Wilkinson, S., Read, A., 2007. Sources of sediment and nutrients to the Gippsland Lakes assessed using catchment modelling and sediment tracers. CSIRO Land Water Sci. Rep. 70/07, pp. 73. Jakeman, A.J., Letcher, R.A., 2003. Integrated assessment and modelling: features, principles and examples for catchment management. Environ. Model. Softw. 18, 491–501. Johnson, I.R., Lodge, G.M., White, R.E., 2003. The sustainable grazing systems pasture model: description, philosophy and application to the SGS national experiment. Aust. J. Exp. Agric. 43, 711–728. Johnson, I.R., Chapman, D.F., Snow, V.O., Eckard, R.J., Parsons, A.J., Lambert, M.G., Cullen, B.R., 2008. DairyMod and EcoMod: biophysical pasture simulation models for Australia and New Zealand. Aust. J. Exp. Agric. 48, 621–631. Keating, B.A., Carberry, P.S., Hammer, G.L., Probert, M.E., Robertson, M.J., Holzworth, D., Huth, N.I., Hargreaves, J.N.G., Meinke, H., Hochman, Z., McLean, G., Verburg, K., Snow, V., Dimes, J.P., Silburn, M., Wang, E., Brown, S., Bristow, K.L., Asseng, S., Chapman, S., McCown, R.L., Freebairn, D.M., Smith, C.J., 2003. An overview of APSIM, a model designed for farming systems simulation. Eur. J. Agron. 18 (3–4), 267–288. Ledgard, S.F., Penno, J.W., Sprosen, M.S., 1999. Nitrogen inputs and losses from clover/grass pastures grazed by dairy cows, as affected by nitrogen fertilizer application. J. Agric. Sci. 132, 215–225. Ledgard, G., 2014. An inventory of nitrogen and phosphorous losses from rural land uses in the southland region. In: Currie, L.D., Christensen, C.L. (Eds.), Nutrient management for the farm, catchment and community. Occasional

T. Thayalakumaran et al. / Agricultural Water Management 178 (2016) 37–51 Report No. 27. Fertilizer and Lime Research Centre. Fertilizer and Lime Research Centre, Massey University, Palmerston North, New Zealand, pp. 16. Littleboy, M., Freebairn, D.M., Silburn, D.M., Woodruff, D.R., Hammer, G.L., 1999. PERFECT Version 3.0. A computer simulation model of productivity erosion runoff functions to evaluate conservation techniques. http://www.apsru.gov. au/apsru/Products/Perfect/PERFECT.HTM. McClymont, D., Freebairn, D.M., Rattray, D., Robinson, B., 2007. Howleaky? Exploring Water Balance and Water Quality Implications of Alternative Land Uses A Computer Program. Dept Natural Resources and Mines, Queensland, Australia. McKenzie, N.J., Jacquier, D.W., Ashton, L.J., Cresswell, H.P., 2000. Estimation of Soil Properties Using the Atlas of Australian Soils. CSIRO Land and Water Technical Report 11/00. http://www.clw.csiro.au/aclep/documents/tr11-00.pdf. Miller, S.A., Landis, A.E., Theis, T.L., 2006. Use of Monte Carlo analysis to characterise nitrogen fluxes in agroecosystems. Environ. Sci. Technol. 40, 2324–2332. Payraudeau, S., van der Werf, H.M.G., Vertes, F., 2007. Analysis of the uncertainty associated with the estimation of nitrogen losses from farming systems. Agric. Syst. 94, 416–430. Rassam, D., Littleboy, M., 2003. Identifying the lateral component of drainage flux in hill slopes. In: Post, D.A. (Ed.), Proceedings of the international Congress on Modelling and Simulation, Townsville, Australia, July 1, 183–188. Refsgaard, J.C., van der Sluijs, J.P., Hojberg, A.L., Vanrolleghem, P.A., 2007. Uncertainty in the environmental modelling processes-a framework and guidance. Environ. Model. Softw. 22, 1543–1556. Rotz, C.A., Corson, M.S., Chianese, D.S., Montes, F., Hafner, S.D., Coiner, C.U., 2012. The integrated farm systems model. Reference Manual Version 3.6. Pasture Systems and Watershed Management Research Unit. Agricultural Research Service, United States Department of Agriculture. Smith, A.P., Western, A.W., 2013. Predicting nitrogen dynamics in a dairy farming catchment using systems synthesis modelling. Agric. Syst. 115, 144–154.

51

Smith, A.P., Western, A.W., Hannah, M.C., 2013. Linking water quality trends with land use intensification in dairy farming catchments. J. Hydrol. 476, 1–12. Snow, V.O., Johnson, I.R., Parsons, A.J., 2009. The single heterogeneous paddock approach to modelling the effects of urine patches on production and leaching in grazed pastures. Crop Pasture Sci. 60 (7), 691–696. Stott, K., O’Brien, G., Thayalakumaran, T., Vigiak, O., Rees, D., Shambrook, D., Roberts, A., 2013. Development of representative livestock systems for the Moe River catchment. Future Farming Systems Research Division, Department of Environment and Primary Industries. 1 Spring Street, Melbourne 3000. Thayalakumaran, T., Vigiak, O., Beverly, C., Roberts, A., Stott, K., Robinson, B., Freebairn, D., 2013. Integrating two process-based models for assessing dairy system management impacts on N losses. In: Piantadosi, J., Anderssen, R.S., Boland, J. (Eds.), MODSIM2013, 20th International Congress on Modelling and Simulation. Modelling and Simulation Society of Australia and New Zealand, December, 2506–2512. Vachaud, G., Chen, T., 2002. Sensitivity of computed values of water balance and nitrate leaching to within soil class variability of transport parameters. J. Hydrol. 264, 87–100. van der Keur, P., Hansen, J.R., Hansen, S., Refsgaard, J.C., 2008. Unceratinty in simulation of nitrate leaching at field and catchment scale within the Odense River Basin. Vadose Zone J. 7, 10–21. Vigiak, O., Stott, K., Rees, D., Shambrook, D., Norng, S., Roberts, A., 2012. Nutrient sampling to enable point and catchment scale modelling in the Moe River. Future Farming Systems Research Division, Department of Environment and Primary Industries. 1 Spring Street, Melbourne 3000. Vigiak, O., Thayalakumaran, T., Beverly, C., Roberts, A., Stott, K., 2013. Estimating the impact of grazing industry on catchment nitrogen loads of the Moe River catchment. In: Piantadosi, J., Anderssen, R.S., Boland, J. (Eds.), MODSIM2013, 20th International Congress on Modelling and Simulation. Modelling and Simulation Society of Australia and New Zealand, December, 2506–2512. Vogeler, I., Beukes, P., Burggraaf, V., 2012. Evaluation of mitigation strategies for nitrate leaching on pasture-based dairy systems. Agric. Syst. 115, 21–28.