Modeling tillage and manure application on soil phosphorous loss under climate change

Phosphorus (P) losses from non-point sources into receiving water bodies play a significant role in eutrophication. Given their failure to adequately control eutrophication in the Lake Erie, conservation recommendations for agricultural watersheds should be reconsidered, particularly under climate change. Using the Environmental Policy Integrated Climate model, the potential impacts on crop yield, surface runoff, tile drainage, and relevant dissolved reactive phosphorus (DRP) losses from manure-amended corn-soybean rotation plots in the Lake Erie basin were estimated for six tillage methods with different mixing efficiencies and manure broadcast application. These were investigated under twelve different regional and global future climate simulations. Tillage alone proved to have only a minor impact on mean corn yield (± 2%). Climate change led to large uncertainties under the single tillage treatment. As a result of the combined effects of biogeochemical processes (e.g., supply) and hydrological (e.g., transport), strong negative relationships (R2 = 0.98) were found between tillage mixing efficiency and DRP loss in surface runoff, tile drainage, and total DRP loss. The impacts of combined manure application (broadcast) and tillage on crop yield and flow volume were similar as those of tillage alone. With respect to total DRP losses, the effects of labile P content change outweighed those of surface runoff or tile drainage change (hydrologic). This resulted in a change in total DRP losses ranging from − 60% to + 151%, with being closely correlated with decreasing tillage mixing efficiency (R2 = 0.94) from moldboard to no-till. Therefore, rotational tillage should be considered for DRP loss reduction and energy saving.


Introduction
Given its critical role in the prolonged eutrophication of aquatic ecosystems, phosphorus (P) loss from agricultural lands via surface runoff and tile drainage to receiving water bodies has become a serious water quality issue, especially in view of the added issues under climate change. Non-point source P pollution leading to eutrophication has left the public sceptical regarding improvements in water quality and recent conservation efforts to reduce sediment-bound P loadings (Smith et al. 2018). For example, lake Erie's cyanobacterial bloom of 2011 was the worst on record (Daloglu et al. 2012;Michalak et al. 2013), and that of 2014 left 400,000 people without potable water (Smith et al. 2015b). Soil P enrichment brought on by conservation practices, compounded with climate change, has led many to doubt whether the 40% reduction of P loading targets announced by the Canada-US Great Lakes Water Quality Agreement (GLWQA) and the Canada-Ontario Agreement on Great Lakes Water Quality and Ecosystem Health (COA) are feasible (Jarvie et al. 2017). Accordingly, we should reconsider the unintended consequences of conservation recommendations (e.g., the adoption of no-till) (Smith et al. 2018).
The mobilization of P from agricultural lands is tied to both supply (biogeochemical, i.e., fertilizer type, and application timing, amount and method) and transport (hydrologic, i.e., surface runoff and tile drainage) factors (Plach et al. 2018). Motew et al. (2018) indicated that given the increase in precipitation events expected under climate change scenarios, a heightened supply of manure P in the soil would likely exacerbate water quality impairment. Similarly, Michalak et al. (2013) noted the concurrence of intense rainfall events and the dissolved reactive phosphorus (DRP) loss from agricultural lands led to Lake Erie's harmful algal bloom of 2011. Ockenden et al. (2017) also showed that, averaged across three representative catchments in the UK, increased rainfall volume would be the most important contributing factor to projected (2050s) increases in winter P losses.
Besides the impacts of climate change on regional hydrology, the supply factor is also of serious concern. Broadcasting fertilizer or manure without incorporation can result in P stratification, and a subsequent rise in DRP loss in runoff (Motew et al. 2018). Stratification often comes hand-in-hand with the reduced-or no-till field management protocols widely implemented throughout the Lake Erie basin to reduce particulate P and soil erosion losses, but ends up increasing DRP in runoff . With a growing implementation of conservation tillage, the perception exists that manure or fertilizer P must be applied by broadcast to adhere to the conditions of notill (Smith et al. 2018). Jarvie et al. (2017) indicated that although these conservation practices reduced particulate P transport and soil erosion, they may unintentionally increase DRP loads. Rather than merely accelerating soil P stratification (Bullerjahn et al. 2016), reduced-or no-till crop field management can accelerate the development and/or retain soil macropores. This increases the connectivity between the surface and subsurface via preferential flow (Jarvis 2007), thereby increasing soil surface P transport to tile drainage by bypassing sorption/desorption processes within the soil matrix (Smith et al. 2015a). Conventional tillage, such as moldboard plow, can reduce buildup of P in the top soil layer, and therefore has potential to reduce risk of DRP loss to Lake Erie . It can also change water flow pathways, altering the balance of surface runoff and tile flow, thereby influencing P transport dynamics (King et al. 2015).
The Environmental Policy Integrated Climate (EPIC) model is a process-based model capable of predicting a crop's physiological growth process arising under a particular crop management strategy and its response to (sub-)daily time scale weather variables, such as atmospheric carbon dioxide concentration (CO 2 ), precipitation and temperature (Schauberger et al. 2017). The EPIC model, calibrated to local conditions, has been widely used to investigate future climate change scenarios (Ahmed et al. 2017;Folberth et al. 2016;Lychuk et al. 2017a;Rosenzweig et al. 2014;Schauberger et al. 2017;Srivastava et al. 2015;Sultan and Gaetani 2016;Xiong et al. 2016). Drawing on the P subroutine developed by Jones et al. (1984), satisfactory estimates of DRP loss through surface runoff and subsurface tile drainage from agricultural lands in the Lake Erie region have been achieved with EPIC (Wang et al. 2018a(Wang et al. , 2018b(Wang et al. , 2019. In addition, EPIC has been tested under both conservation (Le et al. 2018) and conventional (Lychuk et al. 2017b) tillage regimes. Moreover, various needs have been raised from the diversity of stakeholders toward climate change impacts, such as farmers' interests in field-level adaptation and mitigation (Challinor et al. 2009;Redfern et al. 2012). However, no comprehensive study has tackled the question of the potential effectiveness of combined tillage and manure broadcast in limiting DRP loss under future climate change scenarios. Accordingly, beyond the evaluation of climate change impacts (Wang et al. 2021), our objective was to understand the effect of manure broadcast under a range of tillage practices on crop yields, surface runoff, subsurface tile drainage, and relevant DRP losses under twelve RCM ? GCM (regional and global climate models) simulations.

Field experiments
Field experiments were conducted on long-term plots at Agriculture and Agri-Food Canada's Hon. Eugene F. Whelan Experimental Farm at Woodslee, ON, Canada. The soil is a Brookston clay loam classified as an Orthic Humic Gleysol (Soil Classification Working Group 1998), with 26.4% silt, 48.2% sand and 25.4% clay. The soil's bulk density (q) was 1.27 Mg m 3 , while the field capacity and permanent wilting points were 30.4% and 12.7% H 2 O, respectively. The cropping system was a corn-soybean [Zea mays L.
-Glycine max (L.) Merr.] rotation. Solid cattle manure with 50 kg P ha -1 was incorporated after application. Prior to planting corn, ammonium nitrate and potassium chloride with 200 kg N ha -1 and 100 kg K ha -1 were applied, respectively. Herbicides were applied to corn (1.4 kg ha -1 of Roundup, 1.4 kg ha -1 of Dual II and 1.0 kg ha -1 of Atrazine) and soybean (1.4 kg ha -1 of Roundup, 1.4 kg ha -1 of Dual II and 0.5 kg ha -1 of Sencor). A Triple K cultivator and packer were used in spring before planting of either corn or soybean. Chisel plow tillage was applied after harvest. Regular subsurface tile drainage was present in the two replicated plots.
Details regarding the layout of the 67.1 m long 9 15.2 m wide (approximately 0.1 ha) research plots can be found in Tan et al. (1993) and Wang et al. (2018a). Surface runoff and subsurface drainage flow volumes were automatically recorded by water flow meters as well as sending analog and digital signals. A multi-channel data logger utilized the analog signals to monitor, measure and store water volume on a continue basis. Water samples were collected on a flow volume basis with each auto-sampler containing 24 1-L bottles (ISCO Model 2900) activated by digital signals to capture the dynamics of P forms and concentrations (Tan and Zhang 2011;Zhang et al. 2015a). A 1-L sample of surface runoff or tile drainage water was collected from each plot after 1000 L (surface runoff) and 2000 L (tile drainage) of flow volume during the growing season, and after 3000 L (surface runoff) and 5000 L (tile drainage) of flow volume during the non-growing season, respectively. Thereafter, water samples were vacuum-filtered through a 0.45-lm Millipore membrane, then analyzed for DRP (Zhang et al. 2015b).

Climate data
The future climate data were obtained from the North American Regional Climate Change Assessment Program (NARCCAP) (Mearns et al. 2007). The NARC-CAP provides high resolution climate change simulation data allowing investigation of the uncertainty of regional scale projections of future climate and thereafter generation of climate change scenarios for impacts assessment research (Mearns et al. 2012;Mearns et al. 2009). Projections were generated by dynamically downscaling global climate model (GCM) output to 50-km using different regional climate models (RCM). The projections were then downscaled to the Whelan experimental farm at Woodslee, Ontario (experimental location) using the delta method . Six weather variables were used in this study: radiation, maximum and minimum temperature, precipitation, wind speed, and relative humidity. To obtain potential ensembled future daily weather data sets, different RCM ? GCM combinations were used to generate twelve different climatic simulations (Tables 2, 3, 4, 5, 6, 7, and 8). Daily weather data from all 12 NARCCAP simulations were used. Historical simulations spanning the period of 1971-1999 (centered around 1985) and future simulations spanning 2038-2070 (centered around 2055) were developed. Monthly mean differences in radiation, maximum and minimum temperature, precipitation, wind speed and relative humidity between the full historical and full future simulations were calculated. As observed weather data for the study site were available from 2008 to 2015 (centered around 2012), only 43 years apart from the 2038-2070 mid-century period simulated in NARCCAP, monthly mean climate projections were linearly adjusted by a factor of 43/70 [(2055-2012)/(2055-1985)] (Haasnoot et al. 2015;Wang et al. 2015). Then, each month's mean projected temperature changes (°C) adjusted by a factor were added to the observed historical data for daily minimum and maximum temperatures. Likewise, each month's change projections (%) adjusted by a factor were used to adjust the corresponding month of the observed dataset for precipitation, solar radiation, relative humidity, and wind speed, to obtain future weather projections. The CO 2 was assumed to increase from 394 ppm (centered around 2012) to 556 ppm in the future-centered around 2055, interpolated between 2050 and 2060 CO 2 from the A2 climate change scenario, as used in NARCCAP. More details can be found in Wang et al. (2021).

EPIC model simulation
EPIC is a process-based field-scale model which simulates the physicochemical processes in soil and water under different agricultural management scenarios (Williams et al. 2015). The model includes nine separate components: weather, hydrology, erosion, nutrients, soil temperature, plant growth, plant environment control, tillage, and economic budgets. The functions of EPIC's tillage component include: mix nutrients and crop residues, and change in bulk density, ridge height and surface roughness (Sharpley and Williams 1990).
Crop yields, surface runoff, subsurface tile drainage, and relevant P loss from 2008-2015 were used to calibrate and validate the model (Wang et al. 2021). A parameter-crack flow coefficient of 0.4 was used as a substitute for preferential flow partitioned into cracks. The soil layer below tile drains was set to a saturated hydraulic conductivity of 0.01 mm h -1 to retain all the water above tile drains. Thus, lateral subsurface flow above the tile drain was treated as tile drainage in the study site (Wang et al. 2018a). The coefficients regulating P flux between labile and active pools and regulating P flux between active and stable pools were set to 0.025. Phosphorus partition between runoff and sediment was set to 60. Set to 1.35, the soluble P runoff exponent, from the modified GLEAMS method, made soluble P runoff concentration a non-linear function of organic P concentration in soil layer 1. The coefficient for upward movement of soluble P by evaporation was set to 1. The PSC (P sorption coefficient) was determined based on an equation for calcareous soils. Further details regarding the calibration and validation of EPIC for the study site can be found in Wang et al. (2018a). The flux between the labile and active mineral pools was a function of soil temperature, labile P, the active mineral P and PSC (Sharpley and Williams, 1990).
The calibrated and validated EPIC model was used to simulate the water cycle, DRP loss, and crop growth with 8-year projected future climate simulations (2051-2058) and 556 ppm CO 2 (twelve baselines for twelve RCM ? GCM simulations) (Wang et al. 2021). Here the modified management practices were used to study the impacts on crop yield, runoff volume, and DRP loss at the study site. Model parameters and agronomic practices (except for tillage and manure broadcast scenarios) remained unchanged between the baseline and the imposed management practices scenarios. Given tillage practices had no significant effect on soybean yield (data not shown), simulated corn yields, surface runoff, subsurface tile drainage, and relevant DRP loss were compared with each corresponding baseline values for each corresponding RCM ? GCM simulation to analyze the impacts of tillage and manure broadcast under climate change. The relationship fitting was represented in two order polynomials.
Management practice scenarios (Table 1)  Manure broadcast only The only modification was to change the manure application method from incorporation to broadcast, other management practices remained the same, including tillage after harvest.
Combined manure broadcast and tillage B-NT, manure broadcast and no-till; B-SCP, manure broadcast and stubble chisel plow with mixing efficiency of 0.15; baseline, manure incorporated and chisel plow with mixing efficiency of 0.3; B-CCP, manure broadcast and coulter chisel plow with mixing efficiency of 0.4; B-SubCP, manure broadcast and subsoil chisel plow with mixing efficiency of 0.45; B-DP, manure broadcast and disk plow with mixing efficiency of 0.85; B-MBP, manure broadcast and moldboard plow with mixing efficiency of 0.95.

Impacts of tillage on corn yield
Averaged across twelve RCM ? GCM simulations, tillage management had subtle impacts on corn yield, ranging from an increase of 2% with stubble chisel plow to a decrease of 2% with moldboard plow compared with the baseline (chisel plow with mixing efficiency of 0.3) (Fig. 1a). A negative correlation existed between tillage mixing efficiency and variability in corn yield (R 2 = 0.63, Fig. 1a); however, the effects of tillage alone generated large uncertainties in the output of the twelve RCM ? GCM simulations. For instance, the impacts of moldboard plow on corn yield ranged from -10% to 1%, fall well out of interquartile range (IQR) (-3% -0), and stubble chisel plow increased corn yield by 0 to 5% with an IQR of 0 -3% (Fig. 1a). This means that the effectiveness of tillage can be either magnified or eliminated by climate change (i.e., RCM ? GCM simulations), which led to the weak relationship of yield with tillage efficiencies. For example, the effects of different tillage practices on corn yield could range from -10% to ? 5% under HRM3_gfdl with the lowest precipitation of 850 mm y -1 and highest water stress of 6.3 days for baseline, whereas the use of a moldboard plow increased water stress by 6.6 days ( Fig. 2), which resulted in the largest decrease (-10%) in corn yield. On the other hand, the use of a stubble chisel plow decreased water stress by 2.9 days, which led to an increase of 5% in corn yield. The tillage mixing efficiency had a positive correlation with changes in water stress (R 2 = 0.70, Fig. 2).

Impacts of tillage on surface runoff and tile drainage
Averaged across twelve RCM ? GCM simulations, tillage mixing efficiency showed a positive correlation (R 2 = 0.71, Fig. 1b) with changes of surface runoff (-24% -9%) compared with the baseline. Significant fluctuation in forecasted surface runoff was generated under the twelve RCM ? GCM simulations. For example, under the moldboard plow treatment the change in surface runoff ranged from -3% to 17% with IQR of 1% -8% (Fig. 1b); comparatively, no-till led to a decrease in surface runoff ranging from -23% to -9% with IQR of -21% to 15% (Fig. 1b). Climate change did not affect the pattern of surface runoff responding to tillage mixing efficiency, except for the scenario that moldboard plow decreased surface runoff under HRM3_hadcm3 (-3%) and MM5I_had-cm3 (-2%). Generally, tillage was applied after harvest, such that surface runoff increased in the nongrowing season. However, increased potential evapotranspiration (PET), and especially ET under the HRM3_hadcm3 and MM5I_hadcm3 climate scenarios, led to lesser growing season surface runoff, which outweighed the increase of surface runoff during the non-growing season. Climate change affected the magnitude of tillage mixing effectiveness: e.g., the broadest range in surface runoff under different tillage efficiencies (-27% to 17%) occurred under the CRCM_ccsm simulation with projected precipitation of 889 mm y -1 ; while the narrowest range (-16% to 4%) occurred under HRM3_hadcm3 simulation with a projected precipitation of 961 mm y -1 . Disk plow had the largest impacts on surface runoff with an average increase of 9% (2% -16%). The main reason was that ET under moldboard plow exceeded that under disk plow, which resulted in lower surface runoff. The greatest reduction in surface runoff occurred (mean, -24%, range, -28% to 16%) under the stubble chisel plow treatment and was greater than that under no-till (mean, -18%, range, -23% to 13%). One possible reason was the difference in crop water uptake due to the greater increase in corn yield under stubble chisel plow. Another reason was the increase in surface runoff with increased tillage mixing efficiency, which was closely connected with the decrease in tile drainage (R 2 = 0.82, data not shown). For example, stubble chisel plow increased tile drainage flow volume by an average of 21% (16% -32%), which exceeded no-till's average increase of 18%. This means more water went to surface runoff during the partitioning between surface runoff and infiltration.
A strong negative relationship (R 2 = 0.9, Fig. 1c) was found between tillage mixing efficiency and tile drainage, which showed a rate of change ranged from -31% to 21%, compared with the baseline. High uncertainties were generated under single tillage across the twelve RCM ? GCM simulations. For example, moldboard plow decreased tile drainage flow volume by 21% to 41% with IQR of 28% -32%, while no-till increased tile drainage flow volume by 14% to 24% with IQR of 15% -21%. Based on its wider range of change for tile drainage than surface runoff flow volume, tillage mixing efficiency had a greater impact on tile drainage than surface runoff. Moldboard plow had most negative impacts on tile drainage, decreasing its volume by an average of 31% (range, -41% to -21%); while stubble chisel plow increased tile drainage flow volume by 21% (range, 10% -32%).

Impacts of tillage on DRP loss
Averaged across twelve RCM ? GCM simulations, the baseline DRP loss in surface runoff, and tile drainage were 5.32 and 5.83 kg P ha -1 , respectively. The effect of tillage on DRP loss in surface runoff ranged from -98% to 127% and followed a close negative relationship (R 2 = 0.98, Fig. 1d). Similarly, a negative relationship existed between tillage mixing efficiency and DRP loss through tile drainage (R 2 = 0.98, Fig. 1e); but the range of changes (-28% -14%) was much narrower. A similar negative trend (R 2 = 0.98) was apparent for total DRP loss, which varied in a range of -61% to 68% (Fig. 1f). The DRP loads in surface runoff varied in the reverse manner as did changes of surface runoff according to tillage mixing efficiency. Although there is a positive linear relationship between DRP loss and surface runoff flow volume in EPIC (Wang et al. 2018a), the decrease of DRP loss in surface runoff did not respond to changes in the magnitude of surface runoff flow volume in this study. The labile P content in the 0.15 m soil layer correlated with the tillage mixing efficiency (R 2 = 0.77, Fig. 3a). Compared to the chisel plow with a mixing efficiency of 0.3, a lowered efficiency decreased labile P content up to 35% (e.g., no-till). In contrast, a high mixing efficiency increased labile P content by up to 26% (e.g., moldboard plow).

Combined impacts of manure broadcast and tillage
The effects of combination of broadcast application of manure with tillage on crop yield and flow volume were similar to those of tillage alone. Under the manure broadcast treatment, tillage mixing efficiency showed a negative relationship (R 2 = 0.93, Fig. 4a) with DRP loss in surface runoff. The DRP loss in surface runoff ranged from -97% to 335% according to the tillage mixing efficiency from moldboard plow (0.95) to no-till. Increasing tillage mixing efficiency diminished the uncertainties brought about by climate change and showed a gradually narrowing IQR (Fig. 4a). Increasing tillage mixing efficiency also had a positive effect on labile P in the 0-0.15 m soil layer (R 2 = 0.80, Fig. 3b). The greater effect of higher tillage mixing efficiency on reducing DRP loss in surface runoff, overcame broadcast manure's effects on increasing DRP loss as moldboard and disk plow reduced DRP loss by an average of 97% and 74%, respectively (Fig. 4a). When combined with manure broadcast, tillage mixing efficiency showed a weaker correlation with DRP loss in tile drainage (R 2 = 0.69, Fig. 4b) than with tillage alone (R 2 = 0.98). The combined effects of broadcast manure and tillage led to a total DRP loss ranging from -60% to ? 151% (Fig. 4c) that were closely correlated with tillage mixing efficiency (R 2 = 0.94). Generally, lowered tillage mixing efficiency increased corn yield due to a lowered water stress following the close relationship (R 2 = 0.98, data not shown) between water stress and corn yield. Stubble chisel plow increased corn yield to a greater extent than notill because the former decreased more water stress by 0.4 days. For some simulations, such as ECP2_had-cm3, HRM3_hadcm3 and MM5I_hadcm3, the baseline showed a negligible simulated water stress, with minor to no impacts of tillage mixing efficiency on corn yield. Thus, the greater predicted water stress under the RCM ? GCM baseline simulation, the greater the impacts of tillage mixing efficiency on corn yield. Therefore, the climate change would impose large uncertainties on the effectiveness of tillage management on corn yield.
Tillage mixing efficiency was correlated oppositely to surface runoff and tile drainage flow volume Tillage mixing efficiency disrupting the soil continuity of preferential networks had a greater impact on tile drainage flow volume than that on surface runoff flow volume, with a wider range of change and closer correlation. These results were attributable to a high tillage mixing efficiency leading to the disruption of the soil continuity of the preferential network (Williams et al. 2016). Therefore, preferential flow path, simulated as cracks here, reduced the contribution of rainfall events to tile discharge due to disruption. The EPIC model simulated changes in bulk density with varying tillage mixing efficiencies (Sharpley and Williams 1990). Ridge height is calculated according to plow depth. Since different tillage practices may have a varying influence on characteristics of the pore network (shape, number, continuity and size distribution) (Kahlon et al. 2013), a greater tillage mixing efficiency would likely have a greater impact on the disruption of preferential flow path (Williams et al. 2016). This speculation was consistent with our results that tile drainage flow volume reduced as tillage mixing efficiency rose. While our simulations showed impacts of tillage on tile drainage lasted through the non-growing season, Williams et al. (2016) found that the effects of tillage on tile drainage faded away within less than three weeks. This diminished effect of tillage practices was contributed by raindrop impacts, soil wetting and dry cycles which quickly re-established continuous macropores (i.e., consolidation, sealing, and crack formation) after disruption by tillage (Messing and Jarvis 1993). The dynamic reestablishment of macropores was not simulated in EPIC due to

Tillage mixing efficiency negatively correlated to DRP loss
Strong negative relationships between tillage mixing efficiency and DRP loss were due to the combined effects of biogeochemical (e.g., supply) and hydrologic (e.g., transport) processes. The lower content of labile P in the 0.15 m soil layer was associated with high DRP loss in surface runoff and vice versa. This means that the soil with high tillage mixing efficiency retained most labile P in the 0.15 m layer; while the soil under reduced or no-till retained most DRP in the top 0.02 m layer, leading to a high DRP concentration in surface runoff (70% vs. 127%, respectively) that overweighed the decrease in surface runoff flow volume. This was consistent with the finding that 80% of leached manure P stayed in the top soil (0.02 m) and loss into surface runoff through desorption (Vadas et al. 2007). Similarly, Ulén et al. (2010) also demonstrated that no-till may increase DRP loss in surface runoff and subsurface tile drainage due to stratification of fertilizer P or organic matter resulted from increased crop residues. This was consistent with our study that the labile P remaining near the surface of soil under reduced or no-till also led to greater DRP losses in tile drainage (11% vs. 14%, respectively). Zhang et al. (2017) also showed that DRP concentrations in tile drainage were significantly greater under non-tillage than those under conventional tillage. They proposed that DRP transported to tile drainage was contributed from P-rich surface soil, especially for notill under which precipitation transported DRP by entering cracks and other types of macro-pores, such as earthworm furrows and root channels. However, since EPIC does not allow for the simulation of P transport through cracks (Wang et al. 2018a), tillage mixing efficiency affecting DRP concentration in cracks cannot be reflected here. Only disruption of cracks by tillage was simulated and reflected in decreasing tile drainage. Since simulated DRP loss in tile drainage was also correlated to changes in tile drain flow (R 2 = 0.97, data not shown), the combined effects of tillage on DRP loss via biogeochemical (e.g., supply) and hydrologic (e.g., transport) processes led to a negative correlation between DRP loss in tile drainage and tillage mixing efficiency.

Combined effects of manure broadcast and tillage
Manure broadcast oppositely affected DRP loss in surface runoff and tile drainage Broadcast application of manure affected DRP losses was due to the changes of labile P content in soil. Decreased labile P (mean, -8%) in the 0.15 m soil layer was the main factor leading to DRP loss decreases in tile drainage. With no incorporation, broadcast manure makes P left on the soil surface, which results in a high potential for DRP loss through surface runoff. Michalak et al. (2013) also pointed out that a broadcast application of fertilizer was one of the three management practices, along with fall fertilizer application and conservation tillage, contributing most to DRP losses in runoff. Therefore, scenarios were conducted to study the combined impacts of broadcast application of manure and tillage mixing efficiency on DRP losses.

Combined effects on DRP loss dominated by tillage mixing efficiency
The effects of labile P changes (supply) outweighed changes in the magnitude of surface runoff or tile drainage (hydrologic) flow volume, which dominated in DRP loss, were the main reason under the combined manure broadcast and tillage practices. For example, moldboard plow significantly reducing DRP loss through surface runoff was similar to the findings of Zhao et al. (2001) and Ginting et al. (1998). The main reason was that a higher tillage mixing efficiency increased labile P content in the 0.15 m soil layer by as much as 27% (e.g., moldboard plow, Fig. 3B), rather than leaving labile P on the surface of soil as occurred under broadcast application of manure. Reduced-or no-till showed a synergistic interaction with manure broadcast, which accelerated DRP loss in surface runoff. This concurred with the work of Michalak et al. (2013) who showed that a combination of manure broadcast and conservation tillage increased runoff DRP loss. A decrease of up to 51% in labile P in the 0.15 m soil layer under no-till showed that the more labile P remained near the soil surface could easily be lost in surface runoff. However, there is a concern that this DRP loss reduction could be over predicted under combination of manure broadcast and tillage as EPIC did not consider direct DRP loss from manure and through cracks (Wang et al. 2019).
A greater tillage mixing efficiency (moldboard and disk plow) masked the effects of applying manure by broadcast on DRP loss in tile drainage due to its effect on labile P in soil. For example, combined moldboard plow (26% increase alone) and manure broadcast (8% decrease alone) led to a 27% (Fig. 4b) increase of labile P content in the 0.15 m soil layer. As biogeochemical and hydrologic processes regulate P mobilization in agricultural landscapes (Plach et al. 2018), the decrease in tile drainage (31%) in synergy with the increased labile P content (27%) in the 0.15 m soil layer were the main reasons for a decrease in DRP loss through tile drainage. While labile P dominated in its impact over tile drainage volume, DRP loss in tile drainage was not that closely correlated with tile drainage flow volume (R 2 = 0.52). Even if tile drainage increased by as much as 21%, the synergy between manure broadcast and reduced-or no-till resulted in lesser labile P content (up to 51% less, Fig. 3b) in the 0.15 m soil layer, and a concomitant decrease in DRP loss through tile drainage of up to 14%. This also showed how effects of manure and tillage on labile P were dominant in their influence over DRP loss through tile drainage. This concurred with the study of Plach et al. (2018) showing that P supply was more important than hydrology for P loading in tile drainage when soil test P levels were within a reasonable range. For medium tillage mixing efficiencies (e.g., subsoil chisel plow), the combined negative effects of tillage (-17%) and manure broadcast (-9%) also resulted in a decrease in labile P content (-20%) in the 0.15 m soil layer. Thus, a decrease in labile P combined with 10% decrease in tile drainage led to a 14% decrease in DRP losses through tile drainage.
Our study was consistent with the finding that reduced-or no-till in conjunction with broadcasting manure exacerbated DRP loss into the western Lake Erie basin (Jarvie et al. 2017). Climate change impacts diminished as tillage mixing efficiency increased as IQRs ranges for total DRP loss (Fig. 4c) became narrower. The response of total DRP loss to tillage mixing efficiency was similar to that of DRP loss in surface runoff, which implies the latter's change in response to tillage outweigh that of DRP loss in tile drainage or showed with greater variance. This could also be attributable to the limitation of EPIC in simulating P transport in cracks (Wang et al. 2018a), thereby limiting DRP losses in tile drainage in response to tillage.

Conclusions
Here we showed climate change impacts on DRP loss can be mitigated by increasing tillage mixing efficiency, such as from no-till to moldboard plow. In general, higher tillage mixing efficiency decreased DRP losses via surface runoff and tile drainage, even when tillage was combined with a broadcast manure application. However, energy consumption will be an issue considering conventional tillage. Therefore, rotational tillage (between conventional and conservation tillage) should be considered to balance DRP loss reduction and energy conservation under climate change.
The dynamic reestablishment of continuous macropores (i.e., consolidation, sealing, and crack formation) initiated by raindrop impacts, and soil wetting and dry cycles after tillage was not simulated in EPIC due to the simplified usage of a constant crack coefficient to substitute for preferential flow. Another limitation of EPIC in simulating P transport in cracks also limited DRP loss through tile drainage in response to tillage. Thus, in the future, we should focus on the reestablishment of macropores and nutrients transfer through cracks to improve the precision of simulation using EPIC.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Tables 2, 3 ,4,5,6,7,and 8;Fig. 5. Simulation names are the name of the RCM (in capital letters) followed by the GCM simulation it downscaled (in lowercase), with acronyms as used in the NARCCAP data archive (Mearns et al., 2007)