Suction cup system‐dependent variable boundary condition: Transient water flow and multicomponent solute transport

Suction cups are widely used in agricultural and environmental research and monitoring under the hypothesis that the sample chemistry represents the soil pore water solute composition around the cup location. The objective of this study was to analyze the sampling procedures that lead to the most representative sample for soil aqueous phase composition when using a falling head suction cup. This was achieved by simulating simultaneously the hydraulic and geochemical response of the suction cup sampled soil solution and its immediate surroundings when evacuated by a system‐dependent variable boundary condition. Different soils, water contents, vacuum applications, and suction cup internal volumes, as well as variable hydraulic conductivities of the ceramic cup, were evaluated, and their effects on the sampling rate and sample chemical composition were reported. Model results showed that potential extracted soil solution volume depends on a combination of internal suction cup volume and vacuum applied, independently from soil type or water content. A linear relationship was defined between the ratio of the extracted sample to suction cup volume and the initial applied vacuum, for all simulations. The pH values and general chemistry of the sampled solution were found to be more similar to those in the soil when a porous cup system is filled until hydraulic equilibrium is reached. Following this, a small volume suction cup system with a high initial applied vacuum, which allows for faster sample collection, could be optimal.


INTRODUCTION
Measuring soil solution chemistry is of great importance in environmental and agricultural research and monitoring. Dissolved solutes in the soil solution, including nutrients and pesticides, migrate through the soil and potentially leach into groundwater and/or reach other sources of water. From an agricultural perspective, dissolved solutes in the soil solution are available to plants, and their concentration is important for fertilization scheduling and nutrient use efficiency studies. However, measuring soil solution chemistry is a challenging task because the solution needs to be separated from the soil solid phases without significantly changing its composition. Suction cups have been used for this purpose since their invention by Briggs and McCall in 1904. In their work Briggs and McCall (1904) presented for the first time a method of sampling a "portion of soil moisture with the dissolved substances which it contains." The main advantages of suction cups include low installation costs, low disturbance of the water flow regime (depending on sampling F I G U R E 1 Diagram of suction cup system. The porous tip is located in the soil, the suction device applies vacuum, and the water flows from the soil through the porous tip and into the sample collection bottle time), wide range of sampling volumes, nondestructiveness, the ability to sample at various depths, and the possibility of repetitive sampling. The main disadvantages or concerns include the limited and uncertain sampling volume and the ability of this sample to represent the chemical variability encountered in heterogeneous soil profiles (Dorrance, Wilson, Everett, & Cullen, 1991;Geibe, Danielsson, van Hees, & Lundström, 2006;Schlotter, Schack-Kirchner, Hildebrand, & von Wilpert, 2012;Weihermüller, Kasteel, Vanderborght, Šimůnek, & Vereecken, 2011), as well as sample alteration by the system via filtering, sorption, and gas exchange (Grossmann & Udluft, 1991;Hansen & Harris, 1975;Weihermüller et al., 2007). In a comment paper to Briggs and McCall (1904), published in the same year, the representability of the chemistry of the sampled soil solution was already questioned (King, 1904).
In this study, the term "suction cup" is used for the entire sampling device in which a porous cup tip (hereafter "porous tip") in the soil is connected to a sample collection container via a pipe and a flexible tube (Figure 1). The porous tip is buried in the soil, is permeable to water but not to air (Grossmann & Udluft, 1991), and is mostly made of an oxide ceramic, but also from stainless steel, nylon, glass, polyethylene, polytetrafluoroethylene, or polyvinyl chloride (PVC) (Dorrance et al., 1991;Weihermüller et al., 2007). In this study, we will focus on the ceramic porous tip type. The most common, inexpensive device is one to which vacuum is applied and then sealed. As water flows into the cup, due to Core Ideas • A water flow and solute transport model was combined with a chemical modeling tool. • Physical and chemical processes taking place when using suction cups were modeled. • A small volume and high initial vacuum least disturbed the carbonate chemistry. • For systems in equilibrium, the pH error depends on the initial vacuum applied. • A linear relationship between sampled volume and suction cup design is presented.
the potential difference, the vacuum in the cup decreases until it reaches hydraulic equilibrium, and water stops entering the cup.
The spatial extent to which suction cups affect the soil domain has been widely studied (Narasimhan & Dreiss, 1986;Tseng, van Genuchten, & Jury, 1995;van der Ploeg & Beese, 1977;Weihermüller, Kasteel, Vanderborght, Pütz, & Vereecken, 2005;Wu, Baker, & Allmaras, 1995). The activity domain, or sampling volume, depend not only on the soil properties but also on the water content or irrigation rate, applied suction, size of the whole suction cup system or of the porous tip itself, and porous tip hydraulic properties (Essert & Hopmans, 1998). The sampling rate, the spatial extent of a suction cup, and/or its effect on water flow in the profile do not depend solely on the suction cup or the soil properties and conditions, but on functional parameters such as internal volume and vacuum applied and their subsequent change with time (Grossmann & Udluft, 1991;Weihermüller et al., 2007).
Modeling water flow into suction cups has been of interest for 40 yr already, with the objective of understanding their activity domain, representability and disturbance to the water flow, influence on the water content around them, origin of the sampled solution in the soil profile, and suction cup disturbance to solute transport (Tseng et al., 1995;van der Ploeg & Beese, 1977;Warrick & Amoozegar-Fard, 1977;Weihermüller, Kasteel, & Vereecken, 2006;Weihermüller et al., 2011;Wu et al., 1995). Numerical studies can be used to illuminate different aspects in detailed heterogeneous spatial and temporal water flow and solute transport fields and to avoid more complex experimental setups.
A crucial aspect in numerical evaluations of the behavior in a cup-soil system is the implementation of the complex timevariable boundary condition at the porous tip. Narasimhan and Dreiss (1986) presented a numerical model that combines Richards' equation in a finite element mesh for calculating transient flow of water to a suction cup that has been evacuated and then sealed and Boyle's law for calculating the increase in pressure inside the suction cup as water enters the vacuum space. Wu et al. (1995) fitted an exponential function to the measured suction in the cup as a function of time, and Weihermüller et al. (2011) took a −30-cm offset from the depthaveraged pressure head in a simulation without the suction cup.
Chemical modification of the water sample can occur when the soil solution, in equilibrium with the soil gas phase, equilibrates with the gas phase in the suction cup. For example, the partial pressure of CO 2 inside the suction cup is lower than in the soil air. When the extracted soil solution equilibrates with the low-CO 2 atmosphere, it degasses, causing an increase in pH and potentially allowing calcite to precipitate. This phenomenon has been widely studied, mainly for the carbonate chemistry as well as for acidic soils (Hendershot & Courchesne, 1991;Kaupenjohann & David, 1996;Suarez, 1987;Takkar, Ulrich, & Meiwes, 1987;Zabowski & Sletten, 1991), and it was found that the pH shift of the sampled solution will depend on the ratio of sampled solution to the system's total volume. This ratio is dynamic, as soil solution enters the suction cups, and will depend on physical properties of the system. Therefore, a model that combines both physical and chemical processes is needed to realistically describe common scenarios where suction cups are used (e.g., irrigated agriculture). Different techniques have been suggested at the sampler end to prevent discrepancies between the pH of sampled and in situ soil solution. For example, Kaupenjohann and David (1996) used a syringe system instead of suction cups, and Suarez (1987) suggested a multichambered suction cup system. However, as the main advantages of suction cups are their simplicity and low price, the regular falling head type is widely used in environmental and agricultural experimental studies with the objective of monitoring solute transport.
To study the dynamic behavior of water uptake by the suction cup system and the possible effects on the water composition, it is required to simultaneously simulate the hydraulics and the chemistry in the soil, at the boundary between the suction cup and the soil, and inside the suction cup. The objective of this study is to analyze the sampling procedures that lead to the most representative sample for soil aqueous phase composition when using a suction cup. This study presents, for the first time, a framework to simulate a system-dependent variable boundary condition for a suction cup that has been evacuated and sealed, including multicomponent solute transport and chemical equilibrium of a solution that meets an atmosphere with a different partial pressure of CO 2 in a twodimensional model. Such a model can be used to estimate the sampling rates under different soil and sampling conditions and estimate pH errors depending on the chemistry of the system and the ratio of sampled solution to internal suction cup volume.

Model framework
Water flow and multicomponent solute transport in the soil were simulated using HYDRUS (2D/3D) and the UNSATCHEM module (Šimůnek, van Genuchten, & Šejna, 2016), whereas the chemistry of the sampled soil solution in the suction cup was simulated using PHREEQC (Parkhurst & Appelo, 1999).
The UNSATCHEM model simulates the complex interactions between the soil matrix and the major ions in the soil: Na + , Ca 2+ , Cl -, SO 4 2-, K + , Mg 2+ , and alkalinity. Twodimensional advective-dispersive chemical transport under transient water flow conditions in a partially saturated porous medium is described by the partial differential equation (Šimůnek, Šejna, & van Genuchten, 2012): , and N c is the number of aqueous components. The dispersion coefficient tensor in the liquid phase, D ij , is given by Bear (1972) in Šimůnek, van Genuchten, and Šejna (2011). Cation exchange and selectivity between the aqueous and exchangeable phases is defined by the Gapon selectivity coefficient, under the assumption that the cation exchange capacity is constant and independent of pH (Šimůnek et al., 2012). Calcite and gypsum precipitation and dissolution were considered. An expanded explanation of all chemical species considered in UNSATCHEM simulations can be found in Šimůnek et al. (2012). PHREEQC is a numerical code that solves chemical reactions of a wide range of aqueous solutions interacting with minerals, gases, solutions, exchangers, and sorption surfaces, among other processes. For the example in this study, a solution was defined for the PHREEQC simulations including the same aqueous components as in UNSATCHEM, a fixed temperature of 20 • C, solution volume, and pH. A gas phase was defined by volume and CO 2 partial pressure. The two phases were then equilibrated, and calcite and gypsum were allowed to precipitate.
A MATLAB subroutine was built such that during a given time period (1 min was found to be of a high enough resolution in this study), the boundary condition at the porous tip has a fixed pressure head and water flows into the cup due to pressure difference between the soil and the cup. After each time period of 1 min, the change in vacuum in the porous cup system is calculated using Boyle's law as defined in Equation 2. For isothermal conditions, the absolute pressure multiplied by the volume of air remains a constant. For a time period dt i , were h atm is the absolute atmospheric pressure head (1,035 cm), h i is the initially applied vacuum, V i is the total air volume of the suction cup system, including the sample collection bottle (Figure 1), V i+1 is the volume of air after some water entered the cup, and h i+1 is the remaining vacuum pressure. The change in vacuum (h i+1 ) is then used as a new pressure head at the boundary condition of the suction cup, and a new 1-min water flow simulation is done. The partial pressure of CO 2 inside the suction cup is calculated according to Equation 3 after every time period, as it changes with changing pressure: where pCO 2 is the partial pressure of CO 2 , x i is the mole fraction of CO 2 (which in the current case is atmospheric concentration 400 μmol mol −1 at time zero), and P is the total pressure of the gas mixture. The total pressure of the gas inside the suction cup is calculated as where P atm is the atmospheric pressure, and P v is the pressure head (vacuum) applied in the suction cup. The updated pCO 2 , the cumulative solute concentrations leaving the domain through the suction cup boundary, and the water volume are updated every 10 min in an input file of the PHREEQC model. Therefore, all processes happening at the soil-water-air domain are calculated with HYDRUS (2D/3D)-UNSATCEHM, and the reactions happening inside the suction cup are calculated using PHREEQC. The PHREEQC model calculates equilibrium between the liquid phase and the gas phase, allowing for CO 2 degassing of the solution, and calcite and gypsum precipitation. In every time step, the chemical equilibrium is performed for the cumulative volume of water that entered the cup and the corresponding weighted concentration for each solute. Such an approach allows using the atmospheric x i in every simulation time step and letting PHREEQC calculate any mineral precipitation from the beginning of the sampling event. For each new sampling event, the suction cup is opened and the air inside reaches equilibrium with the atmosphere again.
The actions of running a program, either HYDRUS (2D/3D) or PHREEQC, and writing and reading files were all F I G U R E 2 Two-dimensional axisymmetric finite element mesh domain in HYDRUS (2D/3D) with its boundary conditions. The zoomed area shows the ceramic porous tip (light blue) and three subregions of 1-cm width each done through a MATLAB program (codes available at https: //github.com/iaelraij/Scup-VBC or from the corresponding author).

Simulation domain
A suction cup installed in the middle of a cylindrical container was simulated with an axisymmetrical domain ( Figure 2) in HYDRUS (2D/3D). The defined mesh for the general simulations consisted of 4,050 non-equidistant nodes with higher node density in the vicinity of the suction cup and in the porous tip. Longer preliminary simulations aided to define a sufficient fine mesh to deal with the nonlinear chemistry, sharp concentration fronts, and hydraulic gradients. Solute and water mass balances were recorded and monitored during both the preliminary and presented simulations. Mass balances for the simulation cases were calculated by accounting for all fluxes (irrigation, drainage, and sample solution), as well as initial and final conditions. Values were acceptable with a median and average relative errors of <0.05% for water and solute balances, excepting alkalinity, which had high mass balance error percentage of up to 15% (5% on average) due to very low total alkalinity values.
Two materials were defined: soil (loamy sand, loam, and silty loam) and a porous tip made of a ceramic material (Table 1). Soil hydraulic properties for the loamy sand, as well as transport and exchange parameters for all soils, were taken from Raij, Šimůnek, Ben-Gal, and Lazarovitch (2016). Hydraulic properties for the loam and silty loam were obtained from the soil catalog in HYDRUS. Hydraulic properties for the porous tip were taken from Weihermüller et al. (2006). Porous tip bulk density, solute transport, and exchange parameters were assumed equal to those defined for the soil.  The cation exchange capacity of the porous tip was averaged from the measurements of Grobler, Claassens, and Annandale (2003). The bottom of the porous tip was defined at a depth of 0.25 m, and the total depth of the domain was 0.5 m. Porous tip dimensions were defined as follow: diameter of 2.4 cm, length of 6 cm, and wall width of 0.35 cm. Three subsequent subregions of each 1 cm around the porous tip were defined in the geometry of the simulations for the model results to specify the chemistry around the cup in different spatial resolutions (see details in Figure 1). The lower boundary condition was defined as a seepage face with a pressure of −30 cm, in a circular area with a 4.6-cm radius. The boundary condition at the surface of the porous tip of the suction cup was defined as a variable boundary condition through the subroutine between MATLAB and HYDRUS (2D/3D) coupled with UNSATCHEM. Irrigation water was applied uniformly at the upper boundary and defined as variable flux with a rate of 8 mm d −1 (2 L d −1 ) applied during 1 h at 6:00 a.m. In order to test the system with a range of water qualities used for irrigation, applied water quality alternated from low-salinity water for 2 d to a pulse of brackish water for another 2 d, followed by low-salinity water for the rest of the simulation ( Table 2). The initial conditions for water and dissolved, adsorbed, and precipitated solutes were imported from a 60-d simulation in which only low-salinity water was applied until the drainage Cl − concentrations were constant. The CO 2 soil pore air concentrations were defined to be 10 times the atmospheric concentrations, 4,000 μmol mol −1 .

Parametric cases
Different scenarios defined as parametric variations give an overall picture of the simulation approach and reach some general conclusions related to sampling schemes and chemistry of the sampled soil solution (Table 3). The total volume of the suction cup system and the initial vacuum applied were varied in a three-by-three grid (total volume of the suction cup system = 100, 500, and 5,000 cm 3 and the initial vacuum applied = 250, 500, and 813 cm), resulting in nine simulations for each soil type. The sampled solution was limited by the algorithm to 400 cm 3 because higher sampling volumes are generally undesirable and unrealistic. The geometry of the porous tip was not changed; the different volumes occur in the external part of the suction cup, mainly the sample collection bottle, located aboveground ( Figure 1). All the simulations had the same initial and atmospheric boundary conditions. The behavior of the system under different water content conditions was tested using the loam soil as an example by running three additional simulation cases with half of the applied daily irrigation (1 L d −1 ) and zero irrigation for the suction cup with internal volume of 100 cm 3 and for initial vacuum of −250, −500, and −813 cm. The smallest volume (100 cm 3 ) was chosen for convenience, as lower water contents were expected to require longer simulations in order to achieve hydraulic equilibrium.
The sensitivity of the porous tip saturated hydraulic conductivity (K s ) was tested with three additional simulations in the sandy loam soil. The medium value was based on the value of Weihermüller et al. (2006) (Table 1), a value whereas one order of magnitude lower or higher was taken for the low and high values, respectively. The initial vacuum applied was −500 cm, in order to not exceed the air entry pressure values (Soilmoisture Equipment Corporation, 2000). The total volume of the suction cup device was defined as 500 cm 3 for the varying K S simulations.

Physical processes: Parametric cases
A combination of different initial applied pressures (−250, −500, and −813 cm) and internal suction cup volumes (100, 500, and 5,000 cm 3 ) was simulated, and the dynamics between the sampled soil solution volume and pressure head inside the suction cup as a function of time were analyzed. As water enters the cup, the pressure in the cup increases until hydraulic equilibrium is reached (Figure 3a-3c) and water stops flowing into the cup (Figure 3d-3f). During the first minutes after applying vacuum, the water extraction rates are very similar between the different parametric cases for each soil (Figure 3d-3f). The 5-L suction cup system has a distinctive higher sampling volume than the smaller systems, and sampling stops before equilibrium is reached due to a 400-cm 3 limitation imposed to the simulations. Maximum sampled soil solution volume is distinctive for each initial vacuum-internal volume combination and similar for the different soils. Pressure at hydraulic equilibrium is distinctive for each soil, and varies slightly within simulations.
By imposing lower irrigation amounts for the loam soil, water contents in the soil are lower at the time of sampling. Lower water contents in the soil surrounding the suction cup decreased sampling rates and time to equilibrium but did not largely affect the hydraulic equilibrium final pressure and total sampled volumes (Figure 4).
A linear relationship was found for the ratio of potential sample volume (V s ), defined as the volume of water that can enter the sampling device until hydraulic equilibrium is reached, to the total internal sampler volume (V i ) as a function of the initial vacuum applied, h i ( Figure 5). The simulations performed for V i = 5,000 cm 3 were not included, as they did not reach hydraulic equilibrium. This linear relationship can be derived from Equation 2, where V i+1 is the total volume of air remaining inside the suction cup after the system achieved hydraulic equilibrium, and therefore where V s is the volume of the sampled solution when the system achieves hydraulic equilibrium. From Equations 2 and 5, we therefore derive Equation 6: where h f is the pressure (soil and in the suction cup) at hydraulic equilibrium. Even though h f changes along the different simulations (different soils and water content), this change is small in comparison with h i and h atm absolute values, and therefore one linear relationship with two empirical parameters can be defined for all the simulations with R 2 of .998 (Equation 7, Figure 5): As h i decreases, the h f value becomes more important and Equation 7 becomes less accurate.
This relationship can be therefore used to predict the maximum potential soil solution that can be sampled with a Vadose Zone Journal F I G U R E 3 Pressure inside the suction cup and cumulative volume of sampled soil solution for three different suction cup volumes and three initial applied vacuum values. The TAV stands for "time after vacuum application" F I G U R E 4 (a) Pressure inside the suction cup and (b) volume of sampled solution for the 100-cm3 suction cup at three irrigation levels (normal = 2 L d −1 , low irr = 1 L d −1 , and no irr = 0 L d −1 ) and three initial applied vacuum values (P = −250, −500, and −813 cm). TAV stands for "time after vacuum application" system of a known internal volume and initial vacuum applied, if enough time is allowed. Additionally, these two parameters can be easily modified in order to achieve desired sampled volumes.

Chemical processes: Parametric cases
In addition to the physical-related sampling properties, the coupled HYDRUS, UNSATCHEM, and PHREEQC model F I G U R E 5 The ratio of potential sampling solution volume (V s , cm 3 ) to the total internal volume of the suction cup (V i , cm 3 ) as a function of the vacuum application (h i, cm) for a combination of volumes and pressures. Each point represents one sampling day. Low irr. stands for low irrigation treatment at a rate of 1 L d −1 , and no irr. stands for no irrigation treatment F I G U R E 6 Time to hydraulic equilibrium (Δt) as a function of pH error (ΔpH, difference between soil solution sampled pH and soil pH in a 3-cm radius around the cup) for nine initial pressure (P, cm) and internal suction cup volume (V, cm 3 ) combinations in the sandy loam soil. Each dot represents a daily sampling during a 12-d simulation.
The V5000 examples (internal volume of 5,000 cm 3 ) did not reach hydraulic equilibrium, hence the relatively low Δt, which in this case stands for sampling time and not time to hydraulic equilibrium setup also simulates the chemical equilibrium of the sampled solution with the air phase inside the suction cup, as well as the differences between sampled concentrations and soil resident concentrations. Due to the lower CO 2 partial pressure inside the cup compared with the soil pore air phase, the sampled solution degasses and the pH increases. As already shown by Suarez (1987), the pH error decreases when the F I G U R E 7 ΔpH as a function of initial pressure applied (h i, cm) for the different parametric cases (V, cm 3 stands for internal suction cup volume, LowI stands for low irrigation treatment at a rate of 1 L d -1 , NoI stands for no irrigation treatment, and K s stands for the saturated hydraulic conductivity of the ceramic suction cup tip). Each point represents one sampling day, and the ΔpH was calculated once the system achieved hydraulic equilibrium. All 12 d of simulation are presented for the sandy loam soil, whereas only Day 3 is presented for all other soils and cases. Squares represent systems with internal volume of 100 cm 3 , × represents systems with internal volume of 500 cm 3 . Details of the parametric cases can be found in Table 3 sampled volume to total suction cup volume ratio increases. Following equation 7 it can be therefore hypothesized that the pH error will decrease with lower initial vacuum applied.
The pH difference between the sampled soil solution and the soil adjacent to the cup (ΔpH or pH error) is defined as The extraction domain (actual volume from which soil solution is being sampled) changes with time and according to the soil water content or pressure applied at the suction cup boundary (Weihermüller et al., 2005); therefore, a fixed soil volume in a 3-cm radius around the porous tip was chosen for comparison between the pH in the soil and the pH in the soil solution extracted. The ΔpH vs. the time needed to reach hydraulic equilibrium (Δt) is presented in Figure 6 for nine parametric cases using the sandy loam soil. For systems that were allowed to reach hydraulic equilibrium (V = 100 cm 3 and V = 500 cm 3 ), applied vacuum is the primary driver for ΔpH, independent of how long it took to reach equilibrium (Figures 6 and 7), whereas the total internal volume does not affect ΔpH (empty and solid symbols for internal volumes of 100 and 500 cm 3 , respectively). When the initial applied vacuum is low, ΔpH is larger and other factors such as soil texture play a role in the absolute pH error (Figure 7). For the largest volume (5,000 cm 3 ), hydraulic equilibrium is not  Figure 6). The relationship between applied vacuum and ΔpH changes with solution composition, and it is not the same at every sampling day (Figures 6 and 7).
The soil solution volume and quality depends on a combination of the total internal volume, initial vacuum applied, and extraction time. If the internal volume is large, it is possible to collect a large sample, but it will take longer to achieve hydraulic equilibrium. If hydraulic equilibrium is not achieved, then a larger pH error will be obtained. It has already been advised to take as small of samples as possible, in order to reduce the disturbance of the natural soil water flow (Grossmann & Udluft, 1991). However, results from this study show that sampling a low volume of solution in a large system can cause large differences between the pH values in the sampled solution and in the soil surrounding the cup. The soil volume from which the suction cup extracts water is a transient property that depends on the water content at the porous tip boundary and the volume of the water extracted (Grossmann & Udluft, 1991). This dynamic sampling soil volume, together with soil heterogeneity and sampling time, requires a thoughtful interpretation of the measurement and the representativeness of the water extracted. Weihermüller et al. (2007) makes a similar remark and concludes that this limited representativeness can be improved with state-of-the-art modeling techniques, such as the one developed and applied in our study. Figure 8 shows the pressure head inside the cup and the cumulative volume of sampled soil solution for three simu-lations with different hydraulic conductivities of the porous tip. High and medium K s of the porous tip material have similar sampling rates in the first minutes of sampling, and then the medium K s rate decreases ( Figure 8A). However, the high-K s porous tip sampling rate is higher, and it reaches hydraulic equilibrium earlier than the medium-K s cup. The low-K s porous tip samples soil solution with a different rate than the medium-and high-K s tips, and as expected, a longer extraction time is needed to achieve the same soil solution sampled sampled volume as for the high and medium K s .

Physical processes: Porous tip saturated hydraulic conductivity variations
Porous tip properties, such as K s , affect the sampling rate, whereas the potential sampled volume for a specific soil depends on the combination between the initial vacuum and total volume of the system, as this determines how long a gradient can be sustained. The relationships for sample volume and sampling time, as well as the flexibility of the proposed model framework for changing the soil and porous tip hydraulic properties, soil chemistry, imposed vacuum, and porous tip dimensions, can help to design and plan the sampling system and timing, as well as to interpret the results, by revealing the timeframe in which the solution was sampled.

Chemical processes: Porous tip saturated hydraulic conductivity variations
The pH in the soil solution is affected by the different irrigation water qualities and decreases when the pulse of saline water reaches the porous tip at Day 4 ( Figure 9A). In addition, changes in soil solution chemistry along the 12-d simulation affect the pH error in different ways. Figure 9B presents the differences between the pH of the soil solution in the soil volume of 3-cm radius from the porous tip and the pH of the sampled water after equilibrium with F I G U R E 9 pH values for the soil pore water solution (empty circles) collected with three suction cups (Scup) with different saturated hydraulic conductivity (K s ) porous tip values and pH of the soil solution in a 3-cm radius from the porous tip (solid circles) in the sandy loam soil. All "K s " in the legend refers to the K s of the porous tip. (a) pH values at hydraulic equilibrum for the sampling of each day during 12 d and (b) pH values at every 10-min time step, during one sampling event on Day 4 the air phase inside the suction cup system during a single sampling event, for three systems with different porous tip hydraulic conductivities. The pH of the sampled soil solution increases after equilibration between the sampled soil solution and an air phase with lower pCO 2 . The increase in pH is greater for the suction cup with the lowest tip hydraulic conductivity, due to the lower water flux into the cup ( Figure 8). As soil solution enters the suction cup, the pH decreases, and the difference in pH between the soil and the sampled solution decreases as well. The pH decreases due to the increase in pCO 2 inside the cup when the air phase is compressed by the entrance of soil solution and the larger mass of water involved in the chemical equilibrium process.
This simulation setup also simulated calcite precipitation inside the suction cup, along with changes in the Ca 2+ concentration of the sampled solution. The concentration of the sampled solution changes with time, and the values presented in Figures 10 and 11 represent the solution accumulating inside the cup during the simulation. Changing concentrations during the simulation days, due to the pulse of saline water, caused different sampling patterns, as shown for Days 1, 4, 8, and 12 ( Figures 10 and 11). Calcium concentrations in the sampled solution are lower than the soil solution concentrations due to precipitation of calcite inside the suction cup after degassing ( Figure 10). Such chemical reaction inside the suction cup could potentially clog the porous tip pores and therefore decrease its hydraulic conductivity with time. Differences in sampling rate, mainly for the porous tip with the lowest F I G U R E 10 Calcium concentrations in the sampled soil solution (Scup) and in the soil adjacent to the cup in a 3-cm radius (soil) for four sampling events on Days 1, 4, 8, and 12 in the sandy loam soil (Panels a, b, c, and d respectively). The K s value refers to saturated hydraulic conductivity Vadose Zone Journal F I G U R E 11 Chloride concentrations in the sampled soil solution (Scup) and in the soil adjacent to the cup in a 3-cm radius (soil) for four sampling events on Days 1, 4, 8, and 12 in the sandy loam soil(Panels a, b, c, and d respectively). The K s value refers to saturated hydraulic conductivity. The TAV stands for "time after vacuum application" K s , caused differences in resident soil solution concentration, and therefore in the sampled solution concentrations as well. These differences are not present on the first day of sampling ( Figure 9A) but appear in the following days ( Figure 10B-10D). The concentrations of Ca 2+ in the soil, and of Cl − in Figure 11, were calculated as an average of the concentrations in all the nodes at a 3-cm radius. Future work should address the complexity of the changing sampling volume and a comparison of sampled solution with soil concentrations for the actual sampled volume.
Even though Cl − concentrations inside the suction cup are not affected by any geochemical process (degassing, calcite precipitation), the concentrations in the sampled solution differ from the concentrations measured in the soil solution in the 3-cm radius adjacent to the porous tip ( Figure 11). These differences change during the sampling event and with the varying influx concentrations due to the saline water pulse. This is an example of the dynamic sampling area or volume of the suction cup, as the concentration in the extracted solution is representing different locations in the soil at varying times.
Overall, the results of this study show the effect and magnitude of physical properties of a suction cup system on the chemistry of the sampled solution and its similarity to or difference from the soil solution surrounding the sampling point. The potential use for the developed framework is vast and includes calculating the dynamic suction cup sampling soil volume, triggered sampling according to specific conditions around the cup, testing different sampling times, simulating small differences in tracer concentrations (Cl − ) due to concentration gradients close to the soil cup, and so on.

CONCLUSIONS
A combined HYDRUS (2D/3D)-UNSATCHEM and PHREEQC routine enables to simulate water flow and multicomponent solute transport in a soil-suction cup system with a system-dependent boundary condition at the suction cup tip via a MATLAB subroutine. The model simulates the pressure increase in the suction cup when water enters the sampler at rates that depend on the soil hydraulic properties, soil water content, applied sampling vacuum, suction cup volume, and porous tip properties. Sampled solution pH values change during the sampling event due to the chemical equilibrium between the sampled solution and the gas phase inside the suction cup system. From a physical point of view, a small sample volume collected during the shortest time possible has the smallest disturbance in the soil surroundings of the cup. However, pH values and general chemistry of the sampled solution are closer to those in the soil when a porous cup system is filled to its potential. In this case, according to the simulations, a small volume of the suction cup system with a high initial applied vacuum was found to be optimal. The higher the porous tip K s , the higher the sampling rate. However, the K s will be limited by the air entry value of the porous material and, therefore, the maximum possible