**1. Introduction**

Irrigation often results in return flows to streams, either by drainage systems established to remove excess irrigation water in the root zone (drainage returns) or by excess irrigation water moving past the root zone, recharging groundwater and then discharging to streams (groundwater returns) [1]. Return flows can be significant in terms of quantity and/or quality and can potentially impact downstream users or aquatic and riparian ecology [2]. Quantifying the magnitude, timing, and quality of return flows to streams is required to inform the managemen<sup>t</sup> of these risks [3,4]. Drainage returns can be measured directly, but modelling is required to quantify groundwater returns [1].

Several processes need to be represented in the modelling in order to estimate the timing and magnitude of groundwater returns [5,6]. Surface and rootzone processes control the amount of water that moves past the rootzone (rootzone drainage) [7]. This water must pass through the remainder of the vadose zone before it recharges the groundwater system (irrigation recharge). Saturated zone processes, including groundwater–surface water dynamics, control the transmission of this flux to streams.

The surface and rootzone processes are comparatively well understood, being an area of focus for agronomic specialists in the interests of devising improved irrigation practices. The saturated zone also receives attention, being an area of focus for hydrologists and hydrogeologists in the interests of water resource management. The deeper vadose zone, being harder to measure and falling between scientific disciplines, is not as well understood.

Where the vadose zone is thin and relatively uniform with depth, there is minimal delay between rootzone drainage and irrigation recharge such that it is reasonable to assume their equivalence. However, where the vadose zone is thick or includes lithological units where perched water tables may form, the timing and magnitude of irrigation recharge can be significantly a ffected. A thick vadose zone means that it can take a long period of time (in the order of years or decades) for changes in root zone drainage (e.g., due to irrigation development or the introduction of e fficiency measures) to propagate to the water table [8–10]. If a layer of relatively low permeability (e.g., a clay unit) occurs within the vadose zone, a perched water table may form above this layer. In these circumstances, the rate of downwards flux will be restricted by the vertical hydraulic conductivity of the low permeability unit and the perched water table may grow vertically and laterally, and some water could be returned to the surface (e.g., via evapotranspiration or via drainage systems) [11–15].

The inclusion of the deeper vadose zone within catchment-scale models is rare, particularly under irrigated settings. Existing vadose zone models, which can be applied at large scales, mostly assume that water moves under gravity (e.g., [16,17]), but this assumption is invalid where perched water tables form [18].

The challenge of including the deeper vadose zone within models has resulted in di fferent approaches to modelling groundwater returns: forward and inverse modelling. The forward method (e.g., [19,20]), aims to model the relationship between irrigation recharge and groundwater returns. The inverse method (e.g., [15,21]) estimates irrigation recharge using observed groundwater level responses or spring discharges and applies the derived irrigation recharge rates to estimate groundwater returns. There are limitations to both methods. The forward method will often involve simplifying assumptions such as ignoring any time lags or losses associated with perching and use uncertain parameters. Errors due to simplifying assumptions and uncertainty are directly transferable to the predicted groundwater returns. The estimation of recharge by inverse method, which is mostly used where there is no appropriate vadose zone model, does not provide unique solutions. It is also di fficult to formulate scenarios to model di fferent irrigation practices because on-ground actions and their influence on rootzone drainage are not directly linked to irrigation recharge and groundwater returns within the modelling approach.

Both the forward and inverse methods may lead to biases in estimating groundwater returns from irrigation and create uncertainty when linking on-ground actions (e.g., irrigation e fficiency improvements) to impacts on the river. Biases can result in several disbenefits for the managemen<sup>t</sup> of groundwater returns. If the risks of groundwater returns are overstated or thought to occur more acutely than in reality, control measures might be overdesigned, too expensive, or place an unnecessarily restrictive burden on irrigation development. If the risks are understated or thought to occur more gradually than in reality, control measures may be insu fficient to avoid unacceptable impacts. Furthermore, if the modelling approach fails to appropriately link on-ground actions with their e ffects, there can be an inequitable sharing of the cost burden to manage risks.

Recognising the limitations of current methods, the authors of [22,23] developed a semi-analytical approach to model vadose zone transfers. These transfer functions include the simulation of perching and allow for the time lags inherent in unsaturated flow.

The objective of this study was to make use of the transfer functions developed by [22,23] and develop an integrated modelling method to quantify the groundwater returns from irrigation for the managemen<sup>t</sup> of salinity in the lower Murray River, Australia. The integrated approach links the various elements of hydrological system involved: (1) the surface/rootzone, (2) the deeper vadose zone, and (3) the saturated zone. The integration of all three elements is intended to reduce the risk of biases that are inherent in current methods by allowing the models to make full use of the datasets available, as well as to create a transparent process to simulate the impact of various managemen<sup>t</sup> actions. The method

must be applicable on a catchment-wide scale and be practical to implement; that is, it should not impose an unnecessary burden on existing model workflows and their data requirements.

#### **2. Study Area Description**

#### *2.1. Salinity Management in the Murray River*

The Murray River is Australia's longest river. Its catchment (including that of the Darling River) spans more than 1 million km2. It provides water to more than 1.5 million households, including the city of Adelaide, and supports a substantial irrigation industry, being Australia's primary food bowl. Along the lower reaches of the Murray River, low salinity river water (in the order of 100–200 mg/L) is used to irrigate a range of horticultural crops across several irrigation districts. Irrigation has significantly altered the hydrology of this semi-arid environment, leading to much higher recharge rates under areas of irrigation. Because the groundwater is naturally highly saline (typically 18,000 to 35,000 mg/L), the higher rates of recharge can mobilise saline groundwater to the river and its floodplains, potentially resulting in deleterious impacts to downstream users and the environment. An extensive policy framework and managemen<sup>t</sup> system is in place to manage the salinity risks [24]. The policy framework relies on an accounting system that assigns salinity e ffects (both debits and credits) to specific on-ground actions. Groundwater models, representing di fferent irrigation districts, are used as the principal assessment tools to quantify groundwater returns and the salinity credits and debits resulting from irrigation practices.

#### *2.2. Study Area Location, Climate, and Hydrology*

The study area was the Loxton and Bookpurnong irrigation districts of South Australia (Figure 1). Annual rainfall is 275 mm and average maximum temperatures range from 15.6 ◦C in July to 31.3 ◦C in February. Irrigation occurs on the "highland" areas that flank the Murray River and its floodplain.

**Figure 1.** Location of the Loxton and Bookpurnong irrigation districts (the red marking on the inset map shows the location of the study area map).

The hydrology of the study area is dominated by the Murray River, which consists of a main river channel and a meandering network of anabranches and wetlands in its floodplain. There are no local tributaries and rainfall runo ff is negligible. Flow in the Murray River is highly regulated by dams and storages upstream and by a series of weirs in its lower reaches that were constructed to aid navigation. The weirs have created a regime of near-constant water levels, interspersed by floods that are irregular but large.

#### *2.3. Landscape and Geology*

The highland is composed of a dune-swale system, oriented east–west, with relatively subdued relief. It has been incised by the Murray River and its floodplain, which is about 20–30 m below the areas of irrigation.

The landscape is underlain by the Murray geological basin—a 550 m thick sequence of Cenozoic marine sediments. The Murray River has formed a trench through the upper sections of this sedimentary sequence. The floodplain trench is infilled by alluvial sediments to be level with marine Loxton Sands unit of the regional sedimentary sequence. The Loxton Sands are capped by the Blanchetown Clay, a widespread lacustrine deposit that underlies most, but not all, of the irrigation district. It is of variable depth and thickness, tending to occur from depths of 2–3 m below surface as a 5 m thick layer of laminated clays. It is overlain by the aeolian Woorinen Formation, which forms the dune-swale system at the surface and hosts sandy loam soils.
