Modeling water fluxes through containerized soilless substrates using HYDRUS

Containerized crop production can enhance plant health and ensure environmental sustainability, yet proper management requires improved understanding of water fluxes and storage within soilless substrates. Numerical simulation tools developed to simulate water movement in porous media, such as HYDRUS‐3D, may help to quantify effective hydraulic properties of soilless substrates but have not yet been tested in this capacity. Therefore, this study had three main objectives: (a) to assess the accuracy of HYDRUS‐3D for simulating water flow through peat‐ and bark‐based soilless substrates by comparing measured and modeled drainage and water storage; (b) to determine sensitivity of model outputs to individual hydraulic parameters; and (c) to compare model parameterization using three laboratory characterization methods (instantaneous profile sorption, instantaneous profile desorption, and evaporation) vs. inverse modeling with HYDRUS. The results showed that the modeled water contents and drainage timing and amounts were most sensitive to saturated volumetric water content (θs) and least sensitive to saturated hydraulic conductivity (Ks). With regard to parameterization methods, the inverse modeling approach provided the most accurate water balance simulations for both substrates, followed by the sorption method. These two methods estimated lower peak water contents and greater drainage compared with simulations parameterized using desorption and evaporation measurements. Overall, the study results showed that the Richards equation, as calculated using HYDRUS, can provide accurate simulations of water flux through containers when properly calibrated, though sorption‐derived parameters may suffice when model optimization is impractical.


INTRODUCTION
Many specialty crop producers of floriculture, woody ornamental crops, vegetables, and small fruit and fruit trees are and industrial waste products (Raviv, 2013). Containerized production may therefore offer better environmental stewardship and plant health compared with traditional cultivation methods (Barrett, Alexander, Robinson, & Bragg, 2016;Caron & Rochefort, 2013). However, improper management in containerized production can seriously affect crop productivity and cause issues such as overuse of water resources and nutrient losses to surrounding water bodies. Having growing media with optimized hydraulic properties is an important way to alleviate these concerns, making it important to develop robust techniques for characterizing substrates used in these systems.
Typical container production uses soilless substrates comprised, in part, of organic (e.g., Sphagnum peat, pine bark) and inorganic (e.g., perlite, sand) materials. Components are blended to create porous and well-drained substrates that provide adequate air-filled porosity (AS) to alleviate what is known as the "container effect." This effect comes from limited drainage when gravitational head is insufficient to overcome matric tensions (h), thus holding water within pores. Among other problems, this condition results in inhibition of root growth (Bilderback, 1980) and increased risk of pathogens and growth stunting (Bradford, Hsiao, & Yang, 1982). Using soilless substrates in combination with lower drainage holes on the bottom of the containers helps to alleviate waterlogged conditions but may also require excess water application to reduce risk of water stress (Mathers, Yeager, & Case, 2005).
Soilless substrates have been well studied using conventional tools, such as the porometer (Fonteno & Bilderback, 1993), that measure "static" physical properties (i.e., those that are assumed to be constant in time). These measurements attempt to assess whether substrates have an adequate balance of air and water after irrigation (Bilderback et al., 2013). However, they also overlook the importance of time-variable substrate properties, such as unsaturated hydraulic conductivity and gas diffusivity, that must be considered when engineering substrates that properly balance water retention vs. aeration (Caron, Pepin, & Periard, 2014;Turunen et al., 2019). Researchers are increasingly measuring changes in substrate physical (Altland, Owen, & Gabriel, 2011;Cannavo, Hafdhi, & Michel, 2011) and hydraulic properties (Kerloch & Michel, 2015) during production cycles or even during single irrigation events (Raviv, Wallach, Silber, Medina, & Krasnovsky, 1999). These measurements tend to emphasize changes in the air/water ratio as pores become rearranged, particles break down, and roots occupy void space (Allaire-Leung, Caron, & Parent, 1999), as well as the importance of hysteresis phenomena affecting hydraulic properties of such substrates (Naasz, Michel, & Charpentier, 2005;Wever, van Leeuwen, & van der Meer, 1996). However, much uncertainty remains about how to (a) best quantify hydraulic properties in soilless substrates and (b) use these properties to understand water fluxes and Core Ideas • HYDRUS-3D accurately simulated water movement through a container substrate system. • Sorption measurements in the instantaneous profile method provided the highest model accuracy. • Model outputs were more sensitive to saturated volumetric water content than other parameters.
storage in container systems. These sources of uncertainty hinder the ability of producers to adjust these properties for optimal plant health and environmental sustainability. Computational models offer the capability to simulate, interpret, and extend laboratory experiments, and thereby improve understanding of how hydraulic properties influence substrate water dynamics. The HYDRUS numerical model (Šimůnek & Senja, 2011) is commonly used in soil science, engineering, and hydrogeology. This computational program solves the Richards equation (Richards, 1931) using a finite element integration scheme to quantify water movement and storage under unsaturated and saturated conditions. Although widely applied in mineral soils, this and similar models have only recently been applied to predict water flow in soilless substrates Fields, Owen, Stewart, & Heitman, 2017). For example, one-dimensional water movement was simulated in both peat and bark substrates during transient and steady-state conditions using HYDRUS-1D (Fields et al., 2017;, whereas two-dimensional flow was successfully modeled in perlite with HYDRUS-2D (Wever, Nowak, de Sousa Oliveira, & van Winkel, 2004). Other examples include predicting water and heat fluxes over a prolonged period (>1 yr) in green roof substrates (Charpentier, 2015); quantifying water movement and hysteresis for a peat substrate in ebb and flow production (Anlauf, Rehrmann, & Scharat, 2012); investigating substrate water dynamics using Sphagnum peat moss as a representative medium (Weber, Iden, & Durner, 2017); and estimating solute transport in soilless substrates (Boudreau, Caron, Elrick, Fortin, & Gallichand, 2009). To date, there have been few attempts to model water movement and storage in pine bark substrates, even though pine bark comprises >60% of rooting media used by nurseries in the United States (Altland & Krause, 2012). Given the importance of both peat and pine bark as substrates, it is critical to develop hydrologic modeling approaches that are reliable for soilless media (Milks, Fonteno, & Larson, 1989).
The overall goal of this research was to evaluate the optimization of van Genuchten-Mualem parameters for simulating water movement through containerized soilless substrates. This goal was achieved through three separate objectives: (a) to determine the accuracy of a three-dimensional HYDRUS model for predicting water movement through peat-and pine bark-based soilless substrate under both hydrated and dry conditions; (b) to assess the sensitivity of HYDRUS model outputs (water content and drainage timing and quantity) to individual hydraulic parameters; and (c) to compare model parameterizations using laboratory characterization methods (instantaneous profile sorption, instantaneous profile desorption, and evaporation) vs. inverse modeling. The results are intended to help researchers and producers better characterize water drainage and retention dynamics in soilless substrates.

Substrate preparation
Model analysis and the supporting experiments were based on two soilless substrates: (a) a "stabilized" (aged 4-6 mo) pine bark (Pinus taedea L.) that was passed through a 12.6-mm screen (Pacific Organics) and then blended with concrete sand (Heard Aggregates) at a 9:1 volume ratio; and (b) a Canadian Sphagnum peat moss (Fafard) that was hydrated to a mass wetness of ∼3.0 g g −1 and then allowed to equilibrate for 24 h before being blended with horticultural-grade super coarse perlite (Whittemore Company) at a 3:1 volume ratio. The pine bark medium was used to represent conventional open-air nursery production, whereas the peat substrate was used to represent greenhouse production in the southeastern United States.

Substrate hydraulic properties
Static substrate physical properties, including air space (AS), container capacity (CC), total porosity, and bulk density were measured on three replicates of each substrate using a porometer analysis (Fonteno & Bilderback, 1993;Fonteno & Harden, 2010). Air space and CC are two measurements often used in soilless substrate science and are, respectively, defined as the minimum volumetric air content and maximum volumetric water content of a soilless substrate after gravitational flow (i.e., free drainage) occurs (White & Mastalerz, 1966).
Particle size distribution of each substrate was measured by passing three oven-dried replicates of ∼100 g through a stacked column of sieves while agitating with a Ro-Tap shaker (Rx-29; W.S. Tyler) for 5 min. Sieve aperture sizes were 6.30, 2.00, 0.71, 0.50, 0.25, and 0.11 mm. These static physical property values were contrasted between the two substrates using a t test in JMP Pro (12.0.1, SAS Institute).
Three laboratory methods were used to characterize substrate hydraulic properties. The first two methods were based upon the instantaneous profile (IP) method (Watson, 1966), modified for containerized soilless substrates. The third was an evaporation method with the Hyprop instrument (METER Group). Note that substrates were hand packed to similar bulk densities in all test configurations.

Instantaneous profile method
For the IP measurements, a 0.2-m 3 sample of each substrate was loosened and poured into 3.90-L nursery containers (18-cm height) and then slightly compacted by hand to mimic nursery and greenhouse filling operations. The substrate height was then leveled to the top of the container. Containers were then irrigated with 1 L of distilled water applied to the substrate surface with a watering can. Three replicates were prepared for each substrate. Each container was allowed to drain for 30 min between each irrigation to reach CC (Paquet, Caron, & Banton, 1993). Final substrate height was measured to be 16.1 ± 0.13 cm for the Sphagnum peat mix and 16.7 ± 0.06 cm or the bark-amended substrate. This step remains important, as the water desorption curve and the unsaturated hydraulic conductivity can be affected by disturbance (Paquet et al., 1993). Saturated hydraulic conductivity (K s ) was determined by measuring water flux from the container bottom with water supplied to the upper surface to maintain a constant ponding depth via Mariotte bottles. Saturated hydraulic conductivity was calculated with Darcy's law as where Q [L 3 T −1 ] is the total flux through the total area normal to the flux A t [L 2 ], L is the final height of the substrate [L], d is the depth of ponding on the substrate surface [L], and R f is a flux correction factor that accounts for the flow restriction imposed by drainage hole configuration [-] (Allaire, Caron, & Gillichand, 1994). Calculations were conducted with L set equal to measured substrate height and d = 3 cm. Despite a container height of 18 cm, an extension collar was taped around the top wall of container to allow a ponding depth of free water of 3 cm. The R f was experimentally determined to be 1.57 following the procedure described by Allaire et al., 1994. Flow was measured every 2 min for a period of 10 d.
Once K s was calculated, substrates were transferred to a tension table for simultaneous determination of water content and matric potential during a drainage and a rewetting sequence (Caron, Xu, Bernier, Duchense, & Tardif, 1998). Water content was determined using a CS-620 15-cm-long probe inserted vertically and connected to a CR-23X data logger (Campbell Scientific). The measured relative permittivity values were converted to volumetric water contents (θ) using a calibration equation developed for peat substrates (Paquet et al., 1993). Matric potential was measured using two tensiometers equipped with a PX26 pressure transducer (OMEGA Engineering). Tensiometers were inserted horizontally at depths of 3 and 12 cm from the bottom surface of the container. The containerized substrate was saturated and allowed to drain to CC. The water potential at the center of container or substrate profile was then ±8 hPa. Discrete estimates of K(h) [L T −1 ] data were calculated by inverting the Darcy-Buckingham equation with a gradient provided between the tensiometers (9 cm) and change in water storage during (via CS-620 probe volumetric water content meter) providing an estimate of flux through the system.
To characterize both hysteretic branches of the unsaturated hydraulic conductivity and water desorption curves, matric potential and water content measurements of the substrate were taken during further drainage on the tension table from container capacity down to −110 hPa and subsequent rehydration back to container capacity at −8 hPa. This range of soil water potential was chosen because it covered the full range of plant available water to obtain maximum productivity (Caron et al., 1998). We therefore note that the majority of relevant properties associated with growing parameters for soilless media occur at matric potential greater than −110 hPa, and this method was not extended beyond this point (Caron, Elrick, Michel, & Naasz, 2006). When all conductivity measurements were completed, the dry bulk density and total porosity of soil samples were determined by oven drying the cores at 105 o C to determine their dry mass.
The moisture characteristic data from the IP measures were further used to fit a constrained van Genuchten -Mualem model (van Genuchten, 1980;Equation 2) for both the peatand bark-based substrates to provide model parameters for comparison: where S e represents effective saturation [-], θ [L 3 L −3 ], θ s [L 3 L −3 ], and θ r [L 3 L −3 ] represent the respective actual, saturated, and residual volumetric water contents, and α [L −1 ] and n are empirical coefficients [-] that affect the shape and scale of the S e function. The measured water tension (h) and volumetric water content (θ) data from the IP analyses were fit to Equation 2 using SAS 9.4 (SAS Institute), and parameters θ s , θ r , α, and n were determined. The parameters were then used to predict where K s was determined from the average of the three replicate measurements as described above, and τ is a pore connectivity parameter [-] that was set to 0.5, based on Mualem (1976).

Evaporation method
Substrate moisture characteristics (through desorption only) were measured via the evaporative method with the Hyprop instrument (METER Group) and following the procedures listed in Fields, Owen, Zhang, and Fonteno (2016). Briefly, 250-cm 3 cores were uniformly packed with each substrate to a target bulk density equivalent to that of the previous porometer and IP analyses. The media were then saturated and allowed to dry via evaporation while h was measured via tensiometer and θ was calculated through mass loss. Measurements ranged from substrate matric potentials of ∼2 to ∼675 hPa in the peat-based substrate and ∼2 to ∼575 hPa in the bark-based substrate. Measured values were input into HypropFit (Bezerra-Coelho, Zhuang, Barbosa, Soto, & van Genuchten, 2018) to calculate the hydraulic parameters θ s , θ r , α, n, and K s . However, unlike for the IP analysis, in the evaporative analysis, the S e and K(θ) functions were analyzed simultaneously with HypropFit, via a nonlinear algorithm to minimize the sum of weighted square residuals between measured and modeled h and θ data (Peters & Durner, 2008).
As with the IP measurement fitting, τ was set to 0.5.

Mass balance measurements
Two lysimeters were constructed by attaching two 853-cm 2 square Plexiglas plates to either side of a load cell (LSP-10, Transducer Techniques); the lysimeters were then connected to a CR1000 data logger (Campbell Scientific) that was programmed to record load cell weight every second. A 4.0-L container (identical dimensions to container used for IP analyses) was filled with substrate and placed inside a funnel (16-cm i.d.) within a 18.9-L bucket that had a window cut into the side (Figure 1). The entire assembly was then placed on one of the lysimeters. The container was irrigated at a rate of 3.4 cm 3 s −1 via a pressure compensated spray stake (Netafim, Plum model). A wire was wrapped around a plastic outflow tube affixed to the bottom of the funnel and was used to transfer effluent to a container on the second lysimeter ( Figure 1). Irrigation was initiated by actuating a solenoid valve. Timing between the start of irrigation and the onset of effluent drainage was manually recorded, and the load cell measurements were used to infer the outflow rate and water storage within the container for the duration of the experiment. We accounted for the delay between water leaving the container and reaching the collection pan by applying water via the spray stake to an empty container 10 times. The mean time for water to reach the collection pan after draining from the container was 3.0 ± 0.03 s. Fallow (i.e., without plants) containers of each substrate (n = 3 replicates) were allowed to air dry through evaporation to an average θ of 0.292 (peat) and 0.263 cm 3 cm −3 (pine Vadose Zone Journal F I G U R E 1 Representation of the mass balance system designed to measure water flux through a container. Funnel was fit into bucket lid and provided level fit for a 4.0-L container. Lysimeters measured water flux in the system (storage) and outflow (drainage) bark). Each replicate was in turn positioned on the mass balance unit and water was applied over ∼240 cm 2 of the substrate surface at a rate of 3.0 cm 3 s −1 , based on the spray pattern of the pressure compensated spray stake. After conducting this initial set of analyses, the containers were hand watered to reach maximum hydration and allowed to drain for 1 h, before the measurements were repeated under hydrated initial conditions (θ = 0.506 cm 3 cm −3 in the peat and 0.460 cm 3 cm −3 in the pine bark). After this second set of measurements, each substrate was dried at 105 • C for 72 h to determine substrate dry mass and to calculate θ of the container based on the measured weights during the experiments.

Numerical model
A computer representation of the containerized substrate matrix was created using a three-dimensional layered domain in HYDRUS 2D/3D. The domain was an inverted conical frustum with an upper diameter of 19.5 cm, a lower diameter of 15 cm, and a height of 18 cm, with 20 equidistant horizontal layers. The model simulated water fluxes through four drain holes on the matrix sides (2-cm diam.) and a central drain hole on the matrix bottom (2-cm diam.) using the seepage face boundary condition, with the seepage pressure head set to h = 0 ( Figure 2). The upper surface boundary conditions were split into no-flux boundaries and atmospheric boundaries that mimic the water distribution pattern observed during a test of the pressure-compensated spray stake (Figure 2). The matrix had a uniform substrate material. F I G U R E 2 Finite element mesh for the three-dimensional container model. Boundary conditions are represented by colored nodes (white is no flux, green is atmospheric, and brown are seepage faces). The atmospheric boundary condition represents the area where water infiltrates the bulk substrate system and was measured via in situ spray stake pattern Two initial conditions were used: (a) initial water contents of θ = 0.292 cm 3 cm −3 for the peat and 0.263 cm 3 cm −3 for the pine bark (dry scenario), and (b) initial water contents of θ = 0.506 cm 3 cm −3 for the peat and 0.460 cm 3 cm −3 for the pine bark, as measured in the experimental containers after allowing for drainage from overhead irrigation (hydrated scenario). Time-variable atmospheric boundary conditions were set with no flux for the first 30 s to allow for any initial water redistribution (Fields et al., 2017). The boundary condition then simulated irrigation water application at a rate of 0.013 cm s −1 , which represents the application rate of the spray stake corrected for measured infiltration area. The time of irrigation varied based on corresponding experimental conditions: peat-based hydrated = 676 s, peat-based dry = 1,353 s, bark-based hydrated = 426 s, and bark-based dry = 1,287 s. After these times, the upper boundary again reverted to no flux until the end of the simulation time. The hydrated models were run to simulate a total period of 3,000 s, whereas the dry models simulated a total period of 4,000 s.

Sensitivity analysis
We analyzed the sensitivity of model outputs to three hydraulic parameters (θ s , α, and K s ), noting that these three parameters were shown in previous work to influence hydrologic simulation outputs (Chen, Jiao, & Li, 2016;Kotlar et al., 2019;Stewart, Lee, Shuster, & Darner, 2017). We used initial parameters determined via the evaporation method and then individually varied α and K s by factors of 0.1, 0.2, 0.5, 2, 5, and 10 and θ s by increments of 0.1 in a one-at-a-time analysis (Šimůnek & Senja, 2011). We note that this procedure was similar to that used by Ma et al. (2011) and Sadhukhan et al. (2019) to evaluate sensitivity in the Root Zone Water Quality Model 2. Model response was then quantified via forward HYDRUS simulations using four output metrics: (a) time between initiation of irrigation to first drainage, (b) total drainage volume, (c) maximum θ reached during experimentation (θ max ), and (d) θ of the system at the final time step (θ final ). Cumulative flux out of the seepage face (drainage) and mean container water content (θ) were also plotted over time to provide visual observations of how parameter manipulation modified container-water dynamics.

Model optimization and comparison with measured properties
To optimize hydraulic parameters, we created a set of twodimensional axisymmetric models in HYDRUS. We then used the inverse solution module to optimize θ r , θ s , K s , α, and n, specifically by fitting to the cumulative drainage and domain-averaged water contents measured under the dry ini-tial conditions. The optimized parameters were then verified using the full-scale three-dimensional models. A full description of the optimization process, including two-dimensional and three-dimensional model comparison and assessment of initial conditions, is presented in the supplemental material.
The final portion of this experiment involved comparing the optimized model parameters with those constrained using only measured hydraulic properties. Drainage and θ curves were developed for the optimized hydraulic model for both substrates under the hydrated and dry scenarios. Each scenario was then recomputed with hydraulic parameters based on measurements of desorption and/or sorption from the IP analyses and the desorption data captured by the evaporative measurements. Including the optimized parameters, this provided four sets of hydraulic model parameters for each substrate. The modeled predictions for initial drainage, maximum storage, cumulative drainage, and final θ were compared with measured values using RMSE.

Substrate hydraulic properties
The respective CC values of the peat-based and the bark-based substrates were 0.826 ± 0.001 and 0.545 ± 0.011 cm 3 cm −3 , and the AS values of the substrates were 0.091 ± 0.003 and 0.245 ± 0.008 cm 3 cm −3 ( Table 1). The increased CC of peatbased substrate is most likely a result of a larger proportion of fine-sized particles and subsequent increased micropore volume (Drzal, Fonteno, & Cassel, 1999). The bark-based substrate was coarser than the peat-based substrate, with a greater proportion of >2-mm particles ( Table 1).
The difference between the IP desorption and sorption curves (i.e., hysteresis) was also greater for the peat-based substrate than for the bark-based substrate (Table 2), a phenomenon also observed by Naasz et al. (2005) and Otten, Raats, Challa, and Kabat (1999). Greater hysteresis in the peat is probably due in part to the smaller proportion of macropores in peat compared with the bark materials (Nkonglo & Caron, 1999;Tsuneda, Thormann, & Currah, 2001). Hysteresis can also result from hydrophobicity and differential wettability associated with air trapped in peat hydrocysts (specialized structures in peat that provide internal porosity and do not readily reabsorb water once drained [Parent & Ilniki, 2003], or that can be caused by repeated drying cycles Michel, Qi, & Charpentier, 2013]) in peat soils.
The evaporative method ( Figure 3) covered a wider range of matric water potentials than the IP method (Figure 4), in which measurements stopped at h = −110 hPa. For the low-tension range, the evaporative method and the IP desorptive curves had similar θ. For example, θ of the peat-based substrate at h = −100 hPa, as measured by the evaporative Vadose Zone Journal T A B L E 1 Substrate physical properties and particle texture analysis for a 3:1 Sphagnum peat/perlite substrate (v/v) and method, was 0.36 cm 3 cm −3 , whereas the IP desorption analysis gave a value of 0.37 cm 3 cm −3 . The θ values for the bark-based substrate at h = −60 hPa were 0.33 cm 3 cm −3 for the evaporative measurement and 0.29 cm 3 cm −3 for the IP desorption analysis. Thus, the evaporative and IP desorption methods appeared to provide subjectively similar results when used to determine the near-saturated portion of the main drainage curve. Some hysteretic effects were also observed in K(h) of the two substrates (Figure 4). Although both the sorption and desorption data of the peat-and bark-based substrate exhibited an exponential relationship, the sorption and desorption K(h) curves in the peat-based substrate were more visually similar and demonstrated less hysteresis than those in the bark-based substrate (Figure 4). This effect may have been caused by the peat-based substrate possessing increased tortuosity and more closed pores, which created inaccessible pore space (Caron & Nkongolo, 2004).

of 14
FIELDS ET AL.

F I G U R E 4 Substrate moisture tension and hydraulic conductivity
(K) or volumetric water content (VWC) data for a 3:1 (v/v) peat/perlite substrate and a 9:1 (v/v) bark/sand substrate. Data were measured using an instantaneous profile method to measure both desorption (drying) and sorption (wetting) data in a 4.0-L container. Data from three replicates are presented in each panel The evaporative and the IP desorption data resulted in similar parameter values (Table 2), which again indicates the similarity of the two desorption methods. However, the model had lower goodness of fit (R 2 of .35-.69 and .57-.77 for peat-and bark-based substrates, respectively) with the IP measurements than the model with the evaporative measurements (R 2 = .99 and .96 for peat-and bark-based substrates, respectively). This result is possibly due to the smaller range of the measures in the IP method.

Sensitivity analysis
Relative sensitivity of the model outputs to three hydraulic parameters (θ s , α, and K s ) was assessed by comparing the response of four metrics used to describe the system (time to initial drainage, total cumulative drainage, θ max , and θ final ). The model outputs showed high sensitivity to θ s ( Figure 5). For example, increasing and decreasing θ s by 0.1 linearly influenced all four metrics. Cumulative drainage for the peat-based substrate changed by ∼50% for every 0.1 shift in θ s and by ∼20% with every 0.1 shift in θ s in the bark-based substrate. In contrast, model outputs had relatively low sensitivity to K s . For example, cumulative drainage and final θ showed little to no change with the magnitude of K s . Increased K s did reduce time to initial drainage and θ max in both peat-and bark-based substrates. The reduction in θ max likely resulted from increases in K s shifting the entire K(h) curve upwards. In other words, K(h) was able to match the constant inflow rate at a relatively lower water content.
The parameter α had little influence on initial velocity, maximum storage, and final θ for the peat-based substrate ( Figure 5). However, increases in α did increase the total drainage for the peat-based substrate. There was also an observed local minimum in the α sensitivity curve for total drainage in the peat-based substrate. Total drainage started to increase when α was reduced to <0.2 times the original value of α. This local minimum was not observed in the bark-based substrate, where total drainage continued to decrease as α decreased. However, θ max , which decreased at 0.2α original , increased relative to the original value at 0.1α original . These variable effects of α on the four metrics suggest that caution should be taken when calibrating α for specific scenarios.

Hydraulic measurement influences
The final optimized models provided the best fits to measured drainage for the initially hydrated peat-based substrate (RMSE = 64.8 cm 3 ) and similar fits as the IP sorption model for the bark-based substrate (RMSE = 236.7 cm 3 for the optimized model vs. 318.6 cm 3 for the IP sorption model, Table 3). The optimized models provided the closest estimates for θ over the simulation period (RMSE = 0.02 cm 3 cm −3 for the initially hydrated peat and 0.01 cm 3 cm −3 for the initially hydrated bark substrate), whereas the IP sorption parameters resulted in RMSE = 0.04 cm 3 cm −3 for both substrates under the initially hydrated conditions. The θ curves for the hydrated initial condition models parameterized using the IP desorption and evaporative measurements were similar to each other, with 0.29 and 0.28 cm 3 cm −3 greater θ max than the observed data for the peat substrate, and slightly increased variation in the bark substrate (0.19 and 0.08 cm 3 cm −3 greater θ max than observed, Figure 6). The IP desorption and evaporative curves of the peat-based and bark-based substrates with initial hydrated conditions gave poor fits to the observed data for all of the measured metrics, with RMSE values >1,085 cm 3 for drainage and>0.23 cm 3 cm −3 for θ (Table 3). The IP sorption method better matched observations for θ and drainage curves compared with the IP desorption. Although initial drainage occurred earlier (75 s earlier in the peat and 88 s earlier in the bark, representing 30 and 90% of the measured time until drainage) in the sorption model compared with observations in the peat substrate (Table 3), the θ max and total drainage Vadose Zone Journal F I G U R E 5 Sensitivity analysis for four model output metrics based on changes in input values for saturated volumetric water content (θ s ), saturated hydraulic conductivity (K s ), and empirical coefficient α.
were within ∼5% of the measured values. A similar θ max and total drainage was measured in the peat substrate between the IP desorption and evaporative model simulations, with the two models most comparable and yet the furthest from the observed data (Table 3). Similar results were observed for the bark substrate, though in this case, the IP desorption and evaporative methods gave different predictions for θ max and cumu-lative drainage volume, with the IP desorption-parameterized model more closely reflecting the observed values (Table 3).
In the initially dry scenarios, the optimized model error increased for the drainage (RMSE = 365.4 cm 3 in the peat substrate and 435.4 cm 3 in the bark substrate, Table 3). However, the models closely matched the observed θ values (RMSE = 0.01 cm 3 cm −3 for the dry peat and 0.02 cm 3 cm −3 T A B L E 3 Metrics used for determining accuracy of hydrologic simulations of peat-and bark-based substrates based on different hydraulic parameterization methodologies. Both substrates were modeled with initial conditions representing maximum water-holding capacity and volumetric water contents at matric tension = −100 hPa. Time until drainage for the optimized scenarios were taken as the time when 1 cm 3 of cumulative drainage occurred for the dry bark). The IP sorption parameters provided the best model fit for drainage from the bark (RMSE = 238.2 cm 3 ), though they did not match the observations from the peat substrate as well (RMSE = 417.9 cm 3 ). The model runs with IP sorption parameters provided relatively close fits to observed θ values (RMSE = 0.06 cm 3 cm −3 for the dry peat and 0.04 cm 3 cm −3 for the dry bark). As in the hydrated scenarios, the IP desorption and evaporative method parameter sets resulted in poor fits to the observed data, with RMSE values all greater than 788 cm 3 for drainage and 0.15 cm 3 cm −3 for θ.
It should be noted that both substrates, when initially dry, began to drain more rapidly than was simulated with any of the model parameterizations ( Figure 6, Table 3). This result may have been caused by uneven wetting fronts and channeling of water (Hoskins, Owen, & Niemiera, 2014;Lafond, Bergeron, Caron, & Rancourt, 2014), which may also explain why the amount of drainage observed during the experiments exceeded drainage simulated by the models. Although it was possible to optimize the parameters to more closely match the observed drainage patterns (by weighting those data more heavily in the optimization scheme), doing so caused the modeled water contents to be much lower than those measured (data not shown). Further, in our procedure, we optimized the parameters using the dry scenario measurements and then tested those parameters against the initially hydrated experiments (for more information about this procedure, please refer to the supplemental material). The same sets of model parameters worked reasonably well with both scenarios, except that the initially hydrated runs better matched observed water content values when θ s Vadose Zone Journal F I G U R E 6 Cumulative drainage and volumetric water contents (VWC, θ) for a 9:1 (v/v) stabilized pine bark/sand substrate and a 3:1 (v/v) Sphagnum peat/perlite substrate, modeled with HYDRUS 3D in a 4.0-L container with initial moisture conditions set to maximum substrate hydration or equal to a pressure head of −100 cm. Curves represent measured data along with results from simulations parameterized using three different laboratory measurements and with inverse modeling. IP, instantaneous profile values were increased by 0.08-0.09 cm 3 cm −3 (Table 2). This finding provides additional evidence that wetting was probably nonuniform within the initially dry substrates.
Altogether, the results from the hydrated and dry scenarios show that in cases where inverse modeling is infeasible or impractical, the IP sorption method may be the best analysis option for parameterizing HYDRUS for use with soilless substrates. The sorption measurements may be less sensitive than the other methods to water held in the immobile water and in "virtual" pores (i.e., pores that are unable to drain at the tension predicted by capillary theory due to the presence of smaller, water-filled pores around it, as described by Hunt, Ewing, & Horton, 2013). Likewise, air that became trapped during the sorption process was likely a main factor in the lower θ s values determined during that measurement. As organic substrates typically have higher proportions of immobile water than mineral soils (Caron, Price, & Rochefort, 2015), the ability to quantify this property is an important result. Moreover, these results suggest that numerical simulations of water flow and storage may be improved by matching the methods used for hydraulic characterization with the main flow process of interest (i.e., for simulations of wetting processes, parameters determined via wetting experiments may be most appropriate), by including explicit information on hysteresis in soil water retention (as in Anlauf, Rehrmann, & Bettin, 2016), or by accounting for both mobile and immobile water domains. HYDRUS offers modules to account for dual porosity and dual permeability; these modules may be more appropriate for describing preferential flow, but also require additional parameterization and thus remain beyond the scope of the current study (Naasz et al., 2005).

CONCLUSIONS
We used HYDRUS 2D/3D to simulate water fluxes and storage within pine bark-sand and peat-perlite soilless substrates. Our analysis showed that HYDRUS can generate accurate simulations of water content changes and drainage through time, but parameterization is a critical step of this process. To this end we evaluated four different methods to determine parameter values. Three of those methods-IP sorption, IP desorption, and evaporation-required laboratory measurements on repacked samples, whereas the fourth-inverse modeling using HYDRUS-fit parameters to observed flow and wetting data during irrigation experiments. Each method provided a distinct set of hydraulic parameters, highlighting differences based on wetting direction (sorption vs. desorption). At the same time, the laboratory measurements maintained equilibrium conditions, whereas the irrigation experiments may have induced nonequilibrium flow. These distinctions suggest that producer decisions, such as the levels of imposed wetting and drying and irrigation rates, may influence the hydraulic characteristics of the substrate. The findings also point to the importance of matching measurements to relevant processes.
As part of our analysis, we also tested the sensitivity of model outputs to variations in three hydraulic parameters: θ s , K s , and α. Model outputs responded more to changes in θ s than to changes in the other two parameters, suggesting that the amount of available pore space is an important parameter in characterizing soilless substrates. However, our analysis also showed that θ s can be a surprisingly difficult parameter to constrain in these growing media, as the different methods led to several distinct θ s values. Specifically, the IP desorption and evaporation methods provided comparable estimates for θ s , as both techniques measured substrate drying from full saturation. However, the IP sorption measurements, in which substrates were rewetted, led to lower estimates of effective pore space. The inverse model parameterizations ended up determining optimal θ s values in between the two extremes and therefore may have better constrained the amount of pore space that contributed mobile water during the infiltration experiments.
This analysis highlights the importance of capturing effective hydraulic properties associated with real-world water flow and storage dynamics, in which substrates undergo repeated wetting and drying cycles. Proper process depiction is particularly important for porous materials that have high proportions of trapped air and immobile water. Of the three tested laboratory methods, the IP sorption best reflected these factors. Therefore, sorption measures of substrate hydraulics may warrant more widespread adoption by people working to characterize soilless growing media.
Finally, using hydrologic models in containerized production systems may allow producers to better schedule irrigation and manage water dynamics (i.e., avoiding excess drainage or inducing water repellency). These models can also serve to analyze potential new substrates prior to conducting growth tests. Substrate assessments can be costly in terms of time and materials, offering the possibility that these techniques can reduce these expenses. In all, numerical models such as HYDRUS may offer a useful means for growers and researchers to work together to better manage water in container systems.