Role of degradation concepts for adsorbing contaminants in context of wastewater irrigation

Marginal water, like treated effluent or other types of wastewater, induces contaminant accumulation in different environmental compartments such as soil and groundwater. With a root zone model and daily rainfall data as input, accumulation in the case of two different degradation concepts was modeled: (I) the solute degrades only in the solution phase, and (II) the solute degrades in both the solution and the adsorption phase. For linear adsorption, both degradation concepts not only yield similar results, but adsorption does not affect the long‐term concentration. This can be readily explained with the governing equations. For nonlinear adsorption, the degradation concept does affect the concentration and fluxes. When the first concept is valid, the long‐term concentration is not affected by adsorption, whereas the short‐term and fluctuations are. For the second degradation concept, the long‐term average concentration and fluxes are also affected by adsorption. A correction factor was derived to correct the analytical solution of the time‐dependent concentration of nonlinearly adsorbing solutes. For all scenarios, reasonable estimates of the long‐term concentration could be obtained with analytical models for most parameter combinations, which can be used to identify contaminant species that should be prioritized for further research.


INTRODUCTION
Wastewater irrigation is practiced to comply with increasing freshwater scarcity due to the growing world population and climate change (UNDP, 2006). Although wastewater irrigation is currently more widespread in arid and semiarid regions (Hamilton et al., 2005;Lawhon & Schwartz, 2006;Shetty, 2004), it has recently also received increasing attention in temperate climates to decrease the pressure on existing freshwater sources (Bixio et al., 2006).
However, reusing wastewater for irrigation also introduces contaminants to the subsurface. A wide range of contaminant species can be found in wastewater, depending on its source and treatment, such as salts (Feigin, Ravina, & Shalhevet, 1991), heavy metals (Bradl, 2004), organic compounds (Walter, Ederer, Först, & Stieglitz, 2000), and emerging contaminants such as pharmaceuticals (Drillia, Stamatelatou, & Lyberatos, 2005;Geissen et al., 2015;Mrozik & Stefańska, 2014). These contaminant species might undergo adsorption and/or degradation to various degrees. Predicting the environmental fate of such contaminants is crucial to assess the long-term impact Vadose Zone Journal of wastewater irrigation on soil, crop, and groundwater quality.
Simplifications concerning degradation and adsorption are common. The debate regarding degradation has concentrated for years on whether or not adsorbed contaminant can be degraded. Ogram, Jessup, Ou, and Rao (1985) showed that this may not be the case for a pesticide, and also regarding other degradable contaminants, the bioavailability or bioaccessibility issue has been debated. Beltman, Boesten, and van der Zee (2008) showed that the impact of adsorption on leaching and on the longterm fate of solutes for single solute applications under steady-flow conditions depends on whether degradation of adsorbed contaminants occurs. Contrary to the scenario considered by Beltman et al. (2008), wastewater irrigation is characterized by multiple contaminant input pulses and unsteady flow conditions, and it is thus uncertain if these conclusions are relevant for the sustainability assessment of wastewater irrigation. Additionally, adsorption is often assumed to be linear with regard to the concentration in solution, whereas in general, it is a nonlinear process and this can induce significant differences with respect to transport (Cirkel, van der Zee, & Meeussen, 2015;van Duijn & van der Zee, 2018).
For contaminants with small retardation factors, significant fluctuations in concentration may occur in response to daily and seasonal variations in rainfall and irrigation. The effect of erratic atmospheric forcing was investigated using a stochastic approach for soil moisture ; Rodríguez-Iturbe, Porporato, Ridolfi, Isham, & Coxi, 1999) and solute transport, such as salt (Shah et al., 2011;Suweis et al., 2010) and nutrients (Botter, Daly, Porporato, Rodriguez-Iturbe, & Rinaldo, 2008;Manzoni & Porporato, 2009). Analytical solutions are feasible when assuming that daily rainfall follows a Poisson distribution, but seasonal variations in precipitation cannot be accommodated then (van der Zee, Shah, & Vervoort, 2014).
In this paper, we investigate the importance of degradation concepts on the environmental fate of contaminants in the context of wastewater irrigation. We therefore distinguish two degradation concepts: one where the solute is degraded only in the solution phase, and one where the solute is degraded in both the solution and the adsorbed phase. This analysis is carried out for both linearly and nonlinearly adsorbing solutes to identify the complexities that may arise because of nonlinearity. We consider erratic weather with seasonal rainfall patterns to determine its importance for the long-term fate of contaminants. Analytical approximations are developed for average water fluxes to address the dominant features of contaminant redistribution.

Core Ideas
• Results differ for degradation in solution alone and in both solution and adsorbed phase. • Long-term concentration is independent of degradation concept for linear adsorption. • Nonlinear adsorption leads to more degradation and lower concentration. • Analytical models generally approximate the long-term concentration within 10%. • Seasonality is important for drainage and cannot be estimated with annual averages.

MATERIALS AND METHODS
We developed a numerical model to calculate the solute concentration under erratic atmospheric forcing. Root zone models have been used to predict soil moisture and solute concentration dynamics (Rodríguez-Iturbe et al., 1999). Such models consider a perfectly mixed zerodimensional root zone and thus only consider the average moisture content and solute concentration. Adhering to previous studies by Vervoort and van der Zee (2008) and Shah et al. (2011), we consider a root zone of thickness r (cm) that may be assumed homogeneous in view of annual plowing. In the case of heterogeneous soil profiles, effective parameters can still be used to accurately model water flow and solute transport (van der Zee & Boesten, 1991;Wildenschild & Jensen, 1999). The groundwater table is located at a fixed depth (cm) below the root zone. The governing equations are the water and solute mass balance equations. Figure 1 shows a schematic representation of the root zone and the water fluxes.

Water balance
The incoming water fluxes in the root zone are the rates of precipitation, irrigation, and capillary rise from the groundwater, respectively. The outgoing water fluxes are the rates of drainage towards the groundwater, evapotranspiration, and runoff. The water mass balance is therefore where ϕ is the porosity (-), is the soil saturation (-), is precipitation (cm d −1 ), is irrigation (cm d −1 ), is capillary rise (cm d −1 ), is drainage to the groundwater (cm d −1 ), is evapotranspiration (cm d −1 ), and o is runoff (cm d −1 ). The soil saturation varies between 0 and 1, denoting completely dry and fully saturated conditions, respectively.
Erratic rainfall events were assumed to be Poisson distributed in previous studies (Rodríguez-Iturbe & Porporato, 2005). Although this allows for analytical expressions for the probability density function of the soil moisture content, seasonal variations in precipitation cannot be accommodated for, despite being important in many climatic regions. We therefore consider daily rainfall records from two locations with different precipitation patterns instead. The first location is Wageningen, the Netherlands (WAG), using data from the Royal Netherlands Meteorological Institute (KNMI, 2019). There are no strong seasonal variations in precipitation throughout the year, but dry periods can occur sporadically ( Figure 2). The second location is Tennant Creek Airport, Australia (TCA). The rainfall records from the Australian Bureau of Meteorology (BOM, 2019) reveal a strong seasonal variation in precipitation ( Figure 2), with a dry season from April to September, and a wet season from October to March. The chosen rainfall records concern a period of 30 yr, spanning from 1 Apr. 1989 to 31 Mar. 2019 for both locations. The annual cumulative rainfall can vary strongly between years. These rainfall records were analyzed to determine if the Poisson assumption was valid for these locations by means of a chisquare test. This was found to be unlikely for both locations, as both p values were close to zero. All precipitation was assumed to be rain. The amount of precipitation that infiltrates into the soil is dependent on the current storage capacity of the root zone. We assume that excess water can leave the system as runoff when the storage capacity of the soil is exceeded, but this was found to be negligible in our case.
All water fluxes other than precipitation are functions of the saturation, using expressions as derived by Vervoort and van der Zee (2008). The water balance has three saturation thresholds. The first threshold is the wilting point w , below which plants cannot take up water. The second threshold is * , below which the transpiration rate of plants is less than the potential transpiration rate. The third threshold is the field capacity f c , above which water drains to the groundwater. Capillary rise is governed by evapotranspiration demand and can therefore not exceed the evapotranspiration flux. If the maximum possible capillary rise exceeds the evapotranspiration flux, we set the actual capillary flux equal to the evapotranspiration flux. The capillary rise flux, as derived by Vervoort and van der Zee (2008), is then given by where 2 = s ∕ϕ r , with s being the saturated hydraulic conductivity of the soil and = α e (ℎ b ∕ ) 2+3∕ being a dimensionless parameter, ℎ b is the bubbling pressure (cm), and α e = 1 + 3∕[2(1 + 3∕ )] (-) and (-) are soil hydraulic shape parameters (Eagleson, 1978). Additionally, 1 = s ∕(ϕ r {1 − exp[β( * − f c )]}), and β is a parameter to fit the hydraulic conductivity function to the exponential model (Rodríguez-Iturbe & Porporato, 2005). Drainage only occurs when the saturation is above the field capacity . The field capacity is defined as f c = ( ∕ℎ b ) −1∕ , and is thus a function of the soil properties and the depth of the groundwater table. According to Teuling and Troch (2005), evapotranspiration can be expressed as where R is the volumetric root fraction (-), β T is a soil moisture stress function (-), T is a light use efficiency parameter (-), LAI is the leaf area index (-), and p is the potential evapotranspiration rate (cm d −1 ), that is, the evapotranspiration rate when water is not a limiting factor. The parameter β T equals zero below the wilting point, as transpiration cannot occur. When the saturation exceeds * , the evapotranspiration is equal to its maximum. In between these two saturation points, evapotranspiration is assumed to be a linear function of the saturation. The expression for β T then yields These fluxes are shown in Figure 3 as a function of the saturation. Irrigation was not included in the work of Vervoort and van der Zee (2008), nor in the work of Shah et al. (2011).
In this model, we included irrigation in a simplistic manner. When the soil water content drops below a specified threshold irr 1 , and no precipitation occurs that day, a certain amount of irrigation water is applied to bring the saturation to the specified value irr 2 . The duration of irrigation is spread out over 1 d. The saturation thresholds for irrigation are defined as Both parameters a (-) and b (-) can vary between 0 and 1. Irrigation starts when the saturation drops below a certain value between the wilting point and the point below which transpiration becomes limited. The saturation is then brought to a point between field capacity and full saturation.

Solute mass balance
We can write a general mass balance for solutes (contaminants) that allows for adsorption and degradation of solutes. Although groundwater can be a source of solutes, we disregard this and assume groundwater capillary rise will only dilute rootzone water. For a negligible contaminant concentration in rainwater and groundwater, the solute mass balance equation is given by where is the solute mass (g cm −2 ), i is the solute concentration in the irrigation water (g cm −3 ), is the solute Vadose Zone Journal concentration in the root zone (g cm −3 ), is the removal by degradation (g cm −2 d −1 ), and is the removal by plant uptake (g cm −2 d −1 ). The total solute mass is the sum of the mass in solution and the mass adsorbed to the solid phase: = r (ϕ + ρ b ), where ρ b is the soil dry bulk density (g cm −3 ), and is the adsorbed mass fraction (-). We assume that adsorption can be described by the Freundlich isotherm, which relates the adsorbed mass to the solute concentration in solution: = F , where F (cm 3n g −n ), and (-) are empirical constants. This reduces to linear adsorption when = 1.
The description of the degradation rate is crucial for this study. There is ongoing debate on the bioavailability of adsorbed contaminants. Bioavailability is generally defined as the amount of solute that is readily available and accessible for uptake by organisms (Riding, Doick, Martin, Jones, & Semple, 2013). Bacteria may be unable to take up solute from the adsorbed phase, and thus adsorption may reduce the bioavailability of the solute. This was found to be the case for some pesticide species (Guo, Jury, Wagenet, & Flury, 2000;Ogram et al., 1985) and pharmaceuticals (Mrozik & Stefańska, 2014), which were degraded only in the solution phase. However, significant degradation in the adsorption phase can occur (Park, Feng, Ji, Voice, & Boyd, 2003). Whether degradation occurs in the adsorbed phase was found to significantly affect the role adsorption plays for the environmental fate of the solute (Beltman et al., 2008). This distinction is often not made, despite significantly affecting the environmental fate of contaminants. We extend the analysis of Beltman et al. (2008) to transient flow conditions, as is characteristic for wastewater irrigation. Therefore, we consider two scenarios: (I) one where the solute degrades only in the solution phase, and (II) one where degradation occurs in both the solution and adsorbed phase with equal degradation rate constants in both phases. when concept I is valid instead, the degradation rate is a linear function of only the mass in the solution phase where μ is the degradation rate constant (d −1 ) under concept I. When concept II is valid, the degradation rate is expressed as a linear function of the total solute mass: with λ the degradation rate constant (d −1 ) under concept II. To compare both degradation concepts, we modify the degradation rate constant such that the total degradation rates in both concepts are equal. For linear adsorption, this yields where ⟨ ⟩ is the saturation averaged over the whole simulation period, and is the average linear retardation factor: Thus, the degradation rate constant under concept II is smaller than in concept I, yet compensated by the additional degradation in the adsorbed phase. For nonlinearly adsorbing solutes, the retardation factor is a function of the concentration. A rate constant scaled by the retardation factor similar to Equation 10 cannot be derived for the nonlinear case, as it would only be valid for a single concentration. Therefore, we use the same expression as for the linear case.
Finally, we require an expression for the plant uptake rate. This is assumed to be a linear function of a rate constant and the evapotranspiration rate, as most of the contaminants under consideration are not nutrients (Dietz & Schnoor, 2001;Mathur, 2004). We assume that plants only remove contaminants from the solution phase. Thus, we obtain = α , where α is the uptake rate constant (-).
Substituting the above expressions for mass, degradation rate, and uptake rate into the solute mass balance equation and rearranging yields where is the concentration-dependent retardation factor (-).

Analytical solutions
If we disregard the complexity of erratic weather and fast fluctuations of the concentration, we can focus on the main features of the mode such as the average solute concentration and fluxes. The solute mass balance equation then becomes analytically tractable. Although it is already known that this yields reasonable results for most field conditions (Kuntz & Grathwohl, 2009), this may not be the case for extreme infiltration events, such as during irrigation or the seasonal rainfall of the TCA climate. We rewrite Equation 12 with average water fluxes and saturation. Thus is the average nonlinear retardation factor, and ⟨⋅⟩ denotes the timeaverage over the entire simulation period of 30 yr. We rewrite Equation 14 using the dimensionless variableŝ= ⟨ ⟩∕(ϕ⟨ ⟩ r ),̂= ⟨ ⟩ ∕(⟨ ⟩ i ), and the Damköhler numbers Da sol = λϕ⟨ ⟩ r ∕⟨ ⟩, Da plant = α⟨ ⟩∕⟨ ⟩, and Da ads = ( 15) The Damköhler numbers can be used to quickly identify the important transport processes, since they indicate the relative importance of the fluxes compared with leaching. Degradation in the solution phase is reflected by Da sol , plant uptake is reflected by Da plant , and degradation in the adsorbed phase is reflected by Da ads . It should be noted the Damköhler numbers only indicate the average importance of the transport processes over the whole simulation period. Strong variations of the actual importance of each flux may occur through time. Equation 15 also suggests that the concentration at steady state is independent of adsorption when there is no degradation in the adsorbed phase (i.e., Da ads = 0). For linearly adsorbing solutes, Equation 15 can be readily solved for the concentration

16) wherê0 is the dimensionless initial concentration and
is the average linear retardation factor, as defined above. Contrarily to the nonlinear retardation factor, the linear retardation factor is not a function of the solute concentration. A general analytical solution for Equation 15 is not available for the case of nonlinear adsorption, except for special cases when = 0.5 and the solute mass in the solution phase is negligible. This is a reasonable approximation for large retardation factors. Equation 14 can then be rewritten as where is the input rate (g cm −2 d −1 ). For a solute which degrades in both the solution phase and the adsorbed phase, ] is a coefficient for the combined effect of leaching and plant uptake and γ = λ is the degradation rate constant. For a solute which degrades only in the solution phase, and γ = 0. Equation 17 is easily integrated (Boekhold & van der Zee, 1991), yielding 1∕2 . Approximating the solute mass as before, we have for the solute concentration When degradation in the solution phase is significant, Equation 19 must be multiplied by a correction factor to obtain the correct concentration at steady state. The derivation of the correction factor is shown in Appendix A and is defined as This correction factor is applied to Equation 19 at every time step, even though it is derived for the steady state situation. All results for the nonlinearly adsorbing solute presented further on are the corrected concentrations.

Model parameters
The parameter values used in the numerical model are presented in Table 1. The soil properties are based on a sandy clay loam soil from Shah et al. (2011), based on standard Australian soils in Neurotheta (Minasny & McBratney, 2002 irrigation threshold for the start of irrigation a = 0.8, meaning the crops are under occasional drought stress, but the evapotranspiration rate is mostly optimal. The threshold b is chosen as 0, because there is already sufficient drainage due to rainfall without having drainage after irrigation events. The minimum leaching fractions are 22 and 44% for TCA and WAG, respectively. The water fluxes and saturation required for the analytical model have been obtained by averaging the numerical results. The initial concentration cannot be set to zero, as the retardation factor is inversely related to the solute concentration. Therefore, we chose an arbitrary small initial concentration. The value for the plant uptake coefficient was chosen such that plant uptake flux is of similar magnitude as the leaching flux. The degradation rate constant varies between 0 and 12.5 yr −1 , corresponding to solute half-lives of infinity to 20 d, such that the wide range of degradation rates found in soils is represented (Stamatelatou et al., 2003;Zhang, Wu, Sun, Ye, & Chen, 2013). The adsorption coefficient ranges from 0 to 1.1781, such that the average linear retardation factor ranges from 1 to 10. For nonlinearly adsorbing solutes, the retardation factor is a function of the solute concentration and is thus highly variable. The TA B L E 2 Average saturation and water fluxes over the whole simulation period of 30 yr, for a groundwater depth of 200 cm and irrigation parameters a = 0.8 and b = 0 solute concentration in the irrigation water is constant at 0.1 g cm −3 . The numerical scheme is given in Appendix B.

Saturation and water fluxes
Before discussing the role of the degradation concept on solute transport, we first consider the water fluxes. The average saturation and water fluxes are shown in Table 2 for both climates. The soil saturation is held above a certain threshold by irrigation. Therefore, the average saturation and several of the water fluxes are similar under both climates. However, their change over time is significantly different. Therefore, we show the water fluxes over time during the final 5 yr of the simulation (Figure 4). At TCA, both the saturation and its variation are much smaller during the dry season than during the wet season. This causes the average saturation to be slightly lower at TCA than at WAG. The dry season also requires frequent irrigation to prevent the saturation to drop to critical levels. The amount of irrigation is the same per event for both climates, but the higher irrigation frequency at TCA results in larger annual irrigation rates than at WAG. The magnitude of the capillary fluxes from the groundwater is limited by the saturation range induced by irrigation. Therefore, the capillary fluxes do not vary strongly between both climates. However, the capillary flux is generally larger at TCA because of the drier conditions. The drainage events are relatively spread out throughout time at WAG, because there is no distinct seasonality in rainfall. At TCA, however, the drainage events are concentrated in the wet season. The drainage rate per event is generally much larger at TCA, since all the rainfall is concentrated in a small timeframe. Surprisingly, the average annual drainage rate is similar for both climates. This means that these two opposing effects cancel each other out in this particular case. The minimum and maximum evapotranspiration rates are equal between both climates, since the saturation is kept within a certain range and equal potential evapotranspiration rates are used. However, the evapotranspiration rate is generally lower at TCA than at WAG, because the crops experience drought stress during the dry seasons. On average, the annual evapotranspiration rates are similar between both climates. The differences in water fluxes of both climates significantly affects the change of the solute concentration in the root zone, which is expanded upon in Supplemental Figure S1.

Nonreactive solute transport
First, we consider the transport of tracers to see how the climate affects solute transport. Thus, we consider a solute that does not adsorb, does not degrade, and is not taken up by plants. The concentration over time is shown in Figure 5 for both climates. Generally, the concentration is higher for the TCA climate. This is no surprise, because the irrigation rate is much higher, whereas the drainage flux is similar to that of WAG. Since irrigation is the only contaminant source in our model, this results in higher concentrations at TCA. The concentration also follows a seasonal pattern at TCA. The concentration rises during the dry season when there is frequent irrigation and the soil water content is lower, and the concentration drops during the wet season, where there is significant drainage. Nevertheless, the minimum and maximum concentrations show strong variation per year, because of interannual differences in rainfall. Under the WAG climate, the concentration is much lower, because not as much irrigation is required. Strong fluctuations of the concentration are observed here due to occasional dry or wet periods, but they do not follow a seasonal pattern.

Role of the degradation concept
As mentioned before, we distinguish two degradation concepts: (I) a solute that only degrades in the solution phase, and (II) a solute that degrades in both the solution phase and the adsorbed phase. We can then define scenarios L-I and II-L for linearly adsorbing solutes, and I-NL and II-NL for nonlinearly adsorbing solutes. Nonadsorbing solutes can then be considered as a special case of either linear and nonlinear adsorption with F = 0. To illustrate the effect of F I G U R E 6 The daily average dimensionless concentration ∕ i over time as calculated by the numerical model (fluctuating lines) and analytical model (smooth lines) for nonadsorbing, linearly adsorbing, and nonlinearly adsorbing solutes. These simulations used a solute half-life of 50 d and adsorption coefficient of 0.5, giving the average linear retardation factor of 5. WAG, Wageningen, the Netherlands; TCA, Tennant Creek Airport, Australia the degradation concept, the concentration over time for a nonadsorbing, linearly adsorbing, and nonlinearly adsorbing solute is shown under both degradation concepts in Figure 6, for select realizations out of the parameter range mentioned above. The probability density functions of the concertation of these solutes are shown in Supplemental Figure S2.
For linearly adsorbing solutes, the concentration over time is nearly identical under both degradation concepts. Although there is additional degradation in the adsorbed phase in concept II, we used a lower degradation rate constant such that the total degradation rate is equal to that in concept I on average. The degradation rates thus may not be equal in both concepts at specific points in time. As a result, there are minor differences of up to 3% between both degradation concepts. The degradation concept thus does not affect the fate of linearly adsorbing solutes as long as the total degradation rates are equal between both concepts. However, the concentration of nonlinearly adsorbing solutes differs strongly between the two degradation concepts. Significantly smaller concentrations occur under concept II, despite us using a lower degradation rate constant to account for the additional degradation in the adsorbed phase. The solute mass is also lower for nonlinear adsorption under concept II (Supplemental Figure S3). However, the reduction of the degradation rate constant according to Equation 10 is based on linearly adsorbing solutes, since an analytical expression for nonlinear adsorption could not be derived. An even smaller degradation rate constant is required to obtain similar results to concept I. This suggests that the retardation factor of the nonlinearly adsorbing solute is larger than that of the linearly adsorbing solute.
Indeed, the retardation factor of nonlinearly adsorbing solutes is larger than that of linearly adsorbing solutes (Figure 7). The linear retardation factor is equal to 5, on average, and only varies slightly because the fluctuating soil water content (Equation 11). Therefore, the linear retardation factor is similar for both climates, which is why only the result of TCA is presented for the linear case. In contrast, the nonlinear retardation factor is also inversely related to the solute concentration (Equation 13), which causes the nonlinear retardation factor to be larger than the linear retardation factor. Additionally, lower solute concentrations lead to larger retardation factors. The concentrations are higher in TCA, since more irrigation is required to prevent the soil from drying out. Therefore, the retardation factors at TCA are lower than at WAG. Likewise, the concentration is lower under degradation concept II, because of the additional degradation in the adsorbed phase. Therefore, the retardation factors are higher under concept II than under concept I. At low concentrations, the retardation factor is also more sensitive to changes in the concentration. This is seen especially under concept II, where the fluctuations of the retardation factor compared with its mean value are larger at WAG than at TCA.
Despite the differences in retardation factor, Figure 6 shows that the long-term average concentration is independent of the adsorption isotherm under degradation concept I. Additionally, the long-term average concentration was also found to be insensitive to the adsorption coefficient (Appendix C). We can see from Equation 14 that when there is no degradation in the adsorbed phase, the retardation factor only occurs in the time-derivative term. This means that at steady state, when the timederivative term equals zero, the solute concentration and fluxes are independent of the retardation factor. However, Equation 14 considers the simplified case with average water fluxes. Our results confirm this also applies to transient water flow conditions. Similar results were presented by Beltman et al. (2008), who found that the leached mass fraction of single solute applications was insensitive to adsorption under steady flow conditions. As we have shown, this also applies to multiple input pulses and transient flow conditions, and not only the leached mass fraction but also the long-term average solute concentration in the root zone is insensitive to the adsorption coefficient and nonlinearity.
Likewise, the long-term solute fluxes are also insensitive to adsorption under degradation concept I. The effect of the degradation concepts and adsorption isotherms on

WAG II-NL
F I G U R E 8 Solute fluxes relative to the contaminant input averaged over the final ten years of the simulation as a function of the degradation rate constant (λ) at Wageningen, the Netherlands (WAG). The results of Tennant Creek Airport, Australia (TCA) are omitted since they are similar. The adsorption coefficient is set such that the linear retardation factor R = 5 the leached solute mass is shown in Supplemental Figure  S4. The magnitude of the solute fluxes during the final 10 yr relative to the solute input is presented in Figure 8 as a function of the degradation rate constant. The results of scenario II-L are omitted, since they are identical to scenario I-L. The plant uptake rate and leaching rate are similar in all scenarios because of the chosen value of the plant uptake coefficient. The fluxes in scenarios I-L and I-NL are very similar, with only minor differences caused by the difference in retardation factor. The fluxes in scenario II-NL are significantly different. The lower concentrations here enhance adsorption and create a positive feedback loop that promotes degradation in the adsorbed phase, whereas the importance of leaching and plant uptake is diminished. For high degradation rate constants, up to 97% of the incoming solute mass is degraded in scenario II-NL, compared with 75-78% in scenarios I-L and I-NL. Despite adsorption having little effect on the long-term average concentration and fluxes, Figure 6 shows that adsorption does affect the time it takes to reach the longterm concentration. Since the annual rainfall varies per year, a true steady state is never reached. Instead, the longterm concentration fluctuates around an average value, which we consider the steady state concentration for this case. The time it takes to reach steady state is shown as a function of the degradation rate constant and the F I G U R E 9 Time it takes to reach the average long-term concentration for the degradation scenarios for both climates. The results of scenario II-L have been omitted, since they are similar to scenario I-L. WAG, Wageningen, the Netherlands; TCA, Tennant Creek Airport, Australia. See the Nomenclature section for description of variables adsorption coefficient for the abovementioned scenarios in Figure 9. For linear adsorption, the adsorption coefficient can be converted to the retardation factor. However, for nonlinear adsorption, this cannot be done, since the retardation factor is also a function of the solute concentration. Generally, it takes more time to reach the equilibrium situation for higher adsorption coefficients. For linear adsorption, equilibrium is reached faster for large degradation rate constants. For nonlinear adsorption, this does not need to be the case, as lower concentrations increase the retardation factor. Since the retardation factor of nonlinearly adsorbing solutes is larger than that of linearly adsorbing solutes, it takes more time to reach steady state for nonlinearly adsorbing solutes. There are no differences between scenarios I-L and II-L, which is why scenario II-L is omitted. For nonlinear adsorption, it generally takes more time to reach steady state in scenario I-NL than in II-NL. Although the retardation factor is larger for II-NL, the equilibrium situation is reached faster due to the larger degradation rate compared with I-NL.
Additionally, adsorption also affects the amplitude of high-frequency fluctuations. The magnitude of the fluctuations with respect to the average concentration may be quantified by the normalized range wherein 90% of the concentration observations fall: ( 95 − 5 )∕⟨ ⟩, with 5 and 95 as the 5th and 95th percentiles, respectively. This concentration range is shown as a function of the degradation rate constant and adsorption coefficient in Figure 10 for the different scenarios. The dashed lines in Figure 10 indicate where the concentration range equals 1, which means that the fluctuations are equal to the mean value. The results of scenario II-L are again omitted due to their similarity to the results of scenario I-L. In general, larger adsorption coefficients reduce the magnitude of the fluctuations due to the buffering effect of adsorption, whereas larger degradation rate constants increase the magnitude of the fluctuations. The concentration range is smaller in scenario I-NL compared with scenario I-L, since the retardation factor of nonlinearly adsorbing solutes is larger. In scenario II-NL, the fluctuations are larger than in the other scenarios. Rapid degradation in combination with small average concentrations result in large normalized concentration ranges. This indicates that strong deviations from the average concentration occur due to the erratic weather conditions, which may affect the environmental risk of wastewater reuse practices. Fluctuations in concentration also affect the variability of solute fluxes, such as the leached solute mass (Supplemental Table S1).  Figure 6 shows that the short-term concentration is approximated well for TCA for all degradation concepts and adsorption isotherms. However, the concentration at early time is underestimated by the analytical model compared with the numerical results. Similar results have been found for the whole range of degradation rate constants and adsorption coefficients. At WAG, the precipitation rate is low during the first 5 yr compared with the 30-yr average annual precipitation. Therefore, more irrigation is required during this period, which leads to a quick rise in the concentration. However, the analytical model uses the 30-yr annual average fluxes. Extreme situations such as these initial dry years are therefore smoothed out. Time series of irrigation (as in this study) that differ significantly in the first years compared with the entire period may therefore bias the prediction accuracy with simple models.

Approximation by analytical solutions
For the long-term concentration, we indicate the agreement between both models by the ratio of the analytical and numerical solutions an ∕ num . The approximation of the long-term concentration is presented as a function of the degradation rate constant and adsorption coefficient in Figure 11. The results of II-L are omitted since they are comparable with I-L. For I-L and II-L, the analytical approximation of the long-term average concentrations lies within 10% of the numerical result for WAG, and 22% for TCA. The poorest approximations are made for nonadsorbing and nondegrading solutes. For such solutes, the concentration fluctuations are larger since they are not attenuated by adsorption. These fluctuations in concentration affect the average concentration as calculated by the numerical model. However, the fluctuations are not included in the analytical model, which causes the discrepancy between the numerical and analytical models. Similar results are seen for I-NL, but the approximations are slightly worse compared with linear adsorption. For II-NL, the error is up to 20% at both WAG and TCA, but within 10% for most parameter combinations. At WAG, the poorest approximations are seen for large degradation rate constants but small adsorption coefficients. Under such conditions, the concentration fluctuations are large (Figure 10), since there is strong degradation and little attenuation by adsorption. However, the concentration cannot drop below zero, and thus the concentration distribution becomes truncated (Supplemental Figure S2). This results in the average concentration being larger in F I G U R E 1 1 Ratio of the analytical and numerical solutions as a function of the degradation rate constant and the adsorption coefficient. The units of F are cm 3 g −1 for linear adsorption and cm 1.5 g −0.5 for nonlinear adsorption. WAG, Wageningen, the Netherlands; TCA, Tennant Creek Airport, Australia. See the Nomenclature section for description of variables the numerical model than in the analytical model, which yields the larger discrepancies between both models. This effect is not seen at TCA, since the higher irrigation rate there yields a larger concentration, and the concentration distribution therefore does not become truncated.
The long-term leaching and plant uptake solute fluxes were also approximated by the analytical model. The approximation of the leached mass fraction is much poorer than that of the concentration (Figure 12). At WAG, the leached mass fraction is predicted well only for nonadsorbing and nondegrading solutes. The error of the approximation is up to 14% for I-L, II-L, and I-NL, whereas for II-NL, the error is up to 35%. The approximation is considerably poorer for TCA for strongly degrading but weakly adsorbing solutes, with errors up to 327%, whereas strongly adsorbing but weakly degrading solutes are still approximated within 10%. Leaching occurs sporadically in the numerical model, resulting in strong peaks during drainage events, whereas there is no leaching otherwise. In the analytical model, the drainage flux is smoothed out over time, resulting in continuous leaching. The leaching peaks are thus averaged out in the analytical model, which causes a large discrepancy between the numerical and analytical results for TCA. The approximation is poorer at TCA because drainage events only occur during the wet season, whereas drainage events are distributed more evenly at WAG. Leaching is therefore strongly seasonal, which is averaged out by the analytical model. Better results may be obtained with more complex analytical models, such as the periodical leaching model of van der Zee, Shah, van Uffelen, Raats, and dal Ferro (2010), but require more detailed knowledge of the water balance. The approximation of the plant uptake mass fraction is similar to that of the solute concentration (Figure 11), and these contour plots are therefore not repeated here. The analytical model thus predicts plant uptake within 10% of the numerical result for most parameter combinations. The plant uptake fraction is approximated well compared with the leached mass fraction, since it is more constant throughout time.
Although the analytical model yields reasonable approximations for most parameter combinations, it does not account for temporal fluctuations in the concentration, despite being significant under certain conditions. Furthermore, the contaminant mass leached to the groundwater was poorly predicted for weakly adsorbing solutes. Thus, not all environmental risks of contaminant species can be identified from the analytical model. Nevertheless, the analytical model proves useful for obtaining estimates of the long-term situation, which can be used to identify F I G U R E 1 2 Approximation of the leached mass fraction during the last 10 yr of the simulation, as a function of the degradation rate constant and the adsorption coefficient. The units of F are cm 3 g −1 for linear adsorption and cm 1.5 g −0.5 for nonlinear adsorption. WAG, Wageningen, the Netherlands; TCA, Tennant Creek Airport, Australia. See the Nomenclature section for description of variables contaminant species that are likely to cause problems with regard to environmental quality standards. Such information can be used to prioritize contaminant species for more detailed investigation, since the list of contaminant species that are present in wastewater is often extensive.

CONCLUSIONS
When reusing marginal water, such as treated effluent, for irrigation, it is important to have guidelines for assessing whether this is also sustainable. Therefore, we investigated this for adsorbing and degrading contaminants. In view of the debate in the literature regarding what proportion of contaminants in soil is biologically available (or accessible) for microorganism-mediated degradation, we explicitly addressed two degradation concepts on the environmental fate of contaminants. These two degradation concepts are (I) that the contaminant degrades only in the solution phase, and (II) that it degrades in both the solution and the adsorption phase. In reality, an intermediate scenario is the most probable. Daily rainfall records were used from the Netherlands and Australia to consider climates with either little or distinct seasonality, as well as some high-frequency variability of rainfall. We made numeri-cal calculations of the water and solute mass balances and compared them with analytical approximations for average water fluxes. When degradation occurs only in the solution phase, the long-term average solute concentration and fluxes are insensitive to both the extent and the (non)linearity of adsorption. However, adsorption still influences the time to reach periodical steady state and the magnitude of the concentration fluctuations. When adsorbed contaminant can also degrade and adsorption is linear, this does not affect the fate of contaminants compared with the first degradation concept. However, this is different for solutes that adsorb nonlinearly, as then both the concentration and the fluxes of nonlinearly adsorbing solutes change significantly. When degradation also occurs in the adsorbed phase, lower solute concentrations are observed due to more contaminant removal by degradation, at the expense of plant uptake and leaching to the groundwater, which are reduced. These results indicate that the degradation concept of nonlinearly adsorbing solutes is an important property that should be determined in experimental research.
The analytical models could approximate the long-term average concentration and solute fluxes within 10% of the numerical result for most parameter combinations. Although these models do not consider complexities such as the concentration fluctuations that might occur due to erratic weather, they may give a sufficiently good idea of which contaminant species are likely to cause problems with regard to environmental quality standards (e.g., depending on their adsorption coefficient and halflife). The large number of contaminants often found in wastewater, as well as limited knowledge regarding their behaviour in the environment, make sustainability analyses of wastewater irrigation challenging. Such analyses, integrated for soil, crop, and groundwater quality, can be done with the relatively simple-to-use approaches advocated here, particularly the analytical model. This is also useful to identify which contaminant species should receive priority for more detailed analyses.

NOMENCLATURE
A mass input rate, g cm −2 d −1 a coefficient for mass input, g cm −3 d −1 B leaching coefficient for nonlinear adsorption, g (n − 1)/n cm −2(n − 1)/n d −1 b hydraulic shape parameter c solute concentration, g cm −3 C an analytical solution concentration at steady state, g cm −3 C i solute concentration in irrigation water, g cm −3 C num numerically determined average concentration, g cm −3 C 0 solute concentration at initial time, g cm −3 C 5 fifth percentile of concentration, g cm −3 C 95 95th percentile of concentration, g cm −3 c T light use efficiency parameter D coefficient for analytical solution with nonlinear adsorption, d −2 Da ads Damköhler number for degradation in adsorbed phase Da plant Damköhler number for plant uptake Da sol Damköhler number for degradation in solution phase E actual evapotranspiration rate, cm d −1 E p potential evapotranspiration rate, cm d −1 f correction factor analytical solution f LAI leaf area index f R volumetric root fraction G parameter for water flow h b bubbling pressure head, cm I irrigation rate, cm d −1 I a coefficient for irrigation management I b coefficient for irrigation management K F Freundlich adsorption coefficient, cm 3n g −n K s saturated hydraulic conductivity, cm d −1 L drainage rate, cm d −1 M total solute mass, g cm −2 M 0 solute mass at initial time, g cm −2 m parameter for water flow m 1 parameter for water flow m 2 parameter for water flow n Freundlich adsorption exponent P precipitation, cm d −1 Q degradation rate, g cm −2 d −1 q adsorbed mass fraction R average retardation factor R c retardation factor R o runoff, cm d −1 s saturation s fc saturation at field capacity s w saturation at wilting point [s 1 irr saturation threshold for start of irrigation s 2 irr saturation threshold for end of irrigation s * saturation below which evapotranspiration is limited T plant uptake rate, g cm −2 d −1 U capillary rise, cm d −1 Z depth to groundwater table, cm Z r root zone depth, cm α plant uptake rate constant α e hydraulic shape parameter β hydraulic shape parameter β T soil moisture stress function γ degradation coefficient for nonlinear adsorption, d −1 θ volumetric soil water content λ degradation rate constant, d −1 ρ b soil dry bulk density, g cm −3 ρ s soil solid phase density, g cm −3 τ time step size, d ϕ porosity

A C K N O W L E D G M E N T S
This research is financed by the Netherlands Organisation for Scientific Research (NWO) under NWO Contract 14299, which is partly funded by the Ministry of Economic Affairs and Climate Policy, and co-financed by the Netherlands Ministry of Infrastructure and Water Management and partners of the Dutch Water Nexus consortium.

C O N F L I C T O F I N T E R E S T
The authors declare no conflict of interest. How to cite this article: Cornelissen P, van der Zee SEATM, Leijnse A. Role of degradation concepts for adsorbing contaminants in context of wastewater irrigation. Vadose Zone J. 2020;e20064. https://doi.org/10.1002/vzj2.20064

APPENDIX A. Correction factor for analytical solution with nonlinear adsorption
The time-dependent solution of the nonlinear solute mass balance equation (Equation 18) is based on the assumption that most of the mass resides in the adsorbed phase. This can lead to incorrect solutions at steady state when this approximation is not valid. Therefore, a correction factor is required to obtain the correct analytical solution. We start with the steady state version of the approximate mass balance equation (Equation 17). This can be written in dimensionless variables as: The difference with Equation A1 is thus the addition of the Da sol term, representing degradation in the solution phase. The analytical model overestimates the solute concentration at steady state when degradation in the solution phase is significant. Thus, we need a correction factor to correct the analytical solution. We solve Equation A3 to The correction factor is obtained by dividing Equation A4 by Equation A2, such that the approximate solution multiplied by the correction factor returns the same solution as Equation A4 at steady state. This yields the expression for the correction factor as stated in Equation 20. This correction factor is applied to the time-dependent solution (Equation 18). For the steady state solution, Equation A4 can be used directly.

B, Numerical scheme
The water mass balance (Equation 1) and solute mass balance (Equation 12) equations are solved numerically using an explicit scheme. The saturation at the new time step is thus calculated based on the water fluxes at the current time step. The discretized water balance equation is then given by A constant time step size of 1 h is used in all simulations.

C. Average long-term concentration
The average concentration during the final ten years is shown in Figure A1. The shape of the contours is similar between WAG and TCA. However, the concentration is much higher at TCA because there is more frequent irrigation. The contours are mostly controlled by the degradation rate constant. For scenarios L-I and I-NL, adsorption has little effect on the long-term average concentration. These results agree with the findings of (Beltman et al., 2008). For II-NL, the adsorption coefficient does affect the concentration because of degradation in the adsorption phase.