Single‐Rate Dual‐Domain Mass Transfer Model: Elucidating Temperature Effects

In this study, the influence of varying temperatures on the transport behavior of conservative solutes in heterogeneous porous media has been investigated. Column flow experiments employing potassium bromide as tracer were conducted at four temperature levels (3°C, 10°C, 20°C, and 30°C) and the measured electrical conductivity (EC) signal was interpreted through inverse modeling. Additional experiments were performed with deuterium‐enriched water at 10°C and 30°C. For those experiments, deuterium isotope ratios were analyzed alongside anion and cation concentrations. Obtained EC‐based breakthrough curves showed measurable differences in both the observed peak values and tailing intensities that could be clearly attributed to variations in the experimental temperatures. The EC‐based results were further corroborated by the measured isotope ratios and corresponding anion/cation concentrations, although measured differences were less pronounced. The model‐based interpretation of the results employed the standard advection‐dispersion equation as well as three alternative variants that were all based on the single‐rate dual‐domain mass transfer (DDMT) approach, but embracing varying coupling intensities between the experiments. For all variants transport parameters were determined for EC, bromide, and deuterium, respectively. The estimated ranges of the transport parameters all point toward a direct correlation between the effective DDMT rates and the experimental temperature. The observed correlation directly follows the Arrhenius relationship, but is weaker than the one describing the temperature dependence of the molecular diffusion coefficients, therefore pointing to a contribution of non‐diffusive components.

model parameters. Many "classic" tracers are typically considered to behave chemically inert in porous media (Davis et al., 1980;Käss, 2004); however, they can still, for example under changing environmental conditions, alter their expected transport behavior (e.g., Käss, 2004). For instance, the transport of solutes in porous media can be modified by the ambient temperature in the subsurface. Here, documented examples include, among others, temperature impacts on the groundwater flow regime in the form of buoyancy effects and variations of the hydraulic conductivity value (e.g., Engström & Nordell, 2016;Harris et al., 2015;. Furthermore, chemical reaction rates (e.g., Brielmann et al., 2009;Prommer & Stuyfzand, 2005;Saripalli et al., 2001), (micro) biological activity (e.g., Derx et al., 2012;Gharabaghi et al., 2015;Stein et al., 2006) as well as intensity of diffusive processes (e.g., Mon et al., 2016) typically correlate with changes in ambient temperature.
Belonging to the latter group of processes, the diffusive mass exchange between mobile and immobile regions of a heterogeneous porous medium affects solute transport characteristics and explains the non-Fickian anomalous transport in such systems (e.g., Feehley et al., 2000;van Genuchten & Wierenga, 1976). The MADE-1/-2 experiments (see, among others, Boggs & Adams, 1992) are well-known field-scale examples that have provided strong evidence for the importance of physical non-equilibrium processes on generating the often-observed tailing effects in solute transport breakthrough behavior. Also, laboratory investigations (e.g., Haws et al., 2007;Knorr et al., 2016;Rolle et al., 2018;Young & Ball, 1998, and references therein) with well-defined porous media systems have demonstrated the relevance of such processes under more controlled conditions and at much smaller scales.
Being a partially diffusive process (van Genuchten & Wierenga, 1976) and, therefore, also bound to the fundamental process of Brownian motion (see Section 2.1), we hypothesize that the mass exchange between mobile and immobile regions may also show a significant and potentially dynamic temperature dependence that could directly impact the transport behavior of both conservative and reactive substances in porous media. However, to date there is, to our knowledge, no experimental evidence for the existence and magnitude of temperature effects, specifically with regard to mass transfer (exchange) rates. If confirmed, we further speculate that temperature effects could infer some currently unconsidered biases in laboratory-obtained transport and reaction parameters, for example as a result of temperature fluctuations in laboratory environments that are not fully temperature-controlled.
To evaluate and quantify the significance of temperature effects, the specific objective of this article is to experimentally determine the magnitude of temperature-related modifications of tracer breakthrough signals. For this, we performed a suite of laboratory-scale column experiments with potassium bromide (KBr) and deuterium oxide ( 2 H 2 O) as tracers at multiple temperature levels. In the experiments, the electrical conductivity (EC) signal was employed for continuous breakthrough monitoring, and supported by discrete sampling for selected experiments. The collected BTC data sets were then further analyzed through inverse numerical modeling exemplarily employing the single-rate dual-domain transport formulation.

Temperature Dependence of Molecular Diffusion
Molecules dissolved in fluids differ from larger suspended particles with respect to their size and the corresponding intensity of the acting forces, while the molecular kinetic theory of heat still applies (Einstein, 1905). The molecular diffusion coefficient D m is inversely proportional to the prevailing temperature and, assuming laminar flow conditions showing low Reynolds numbers, the Stokes-Einstein relation (Einstein, 1905;also, e.g., Miller, 1924) can be employed for calculating this coefficient, that is, whereby k B represents the Boltzmann constant with a value of ∼1.380649 × 10 −23 m 2 kg s −2 K −1 , T is the absolute temperature,  ≈ 3.14159 is Archimedes' constant, r is the mean hydrodynamic radius (here: of the BINDER ET AL.

10.1029/2020WR029474
2 of 20 investigated tracer molecules) and µ is the temperature-dependent dynamic viscosity of the fluid (here: the ambient pore water). In order to consider scenarios with solutes and solvents of similar size, equations have been adopted accordingly (e.g., Glasstone et al., 1941;Sutherland, 1905;as cited in Cussler [2009]), each of them with "T/µ(T)" as central part.   µ T can be calculated empirically by employing Equation 2 (after Voss, 1984), with [T] = K and     µ = kg m −1 s −1 , that is, for waters between 0°C and 30°C.

Concept of Dual-Domain Mass Transfer in Porous Media
Non-Fickian, anomalous transport in natural porous media can be described by a wide range of different approaches (for a comparative overview, see Lu et al., 2018), including the multi-rate mass transfer (MRMT) model. The MRMT model class includes several approaches incorporating non-equilibrium mass exchange. These models are mathematically almost equivalent but differ with respect to their formulation describing the mass transfer process, ranging from the basic single-rate dual-domain model (hereafter referred to as DDMT) to advanced complex multi-domain/multi-rate models including such following power laws (e.g., Babey et al., 2015;Haggerty & Gorelick, 1995;Lu et al., 2018;Pedretti et al., 2013).
For the present study the single-rate DDMT was selected, following the "parsimonious model" approach (see, e.g., De Aguinaga, 2010) in order to avoid over-parameterization of the model. In addition, the DDMT is also the most widely used formulation, not least because it is implemented in several of the most widely used solute transport simulators such as MT3DMS (Zheng & Wang, 1999), MT3D-USGS (Bedekar et al., 2016) and HYDRUS (Šimůnek et al., 2009). The DDMT partitions the sediment's total porosity n into two separate domains: the mobile (n m ) and the immobile domain n im (e.g., van Genuchten & Wierenga, 1976). In the "single-rate" DDMT, the diffusive mass exchange, being proportional to the concentration difference between the two domains and following its gradient, is controlled by a single parameter: the first-order mass transfer rate constant (Coats & Smith, 1964;Valocchi, 1985). Although the single-rate DDMT is clearly a simplification of the complexities impacting real-world aquifers, it has often been utilized to describe the solute transport behavior observed in column experiments as well as in field-scale tracer tests with surprising accuracy.
In the absence of any chemical reactions, the governing transport equations for the mobile (Equation 4) and immobile domains (Equation 5) can be described as follows (1-dimensional (1-D) physical non-equilibrium; e.g., Coats & Smith, 1964;van Genuchten & Wierenga, 1976;and others): and BINDER ET AL.

10.1029/2020WR029474
3 of 20 where C m (t,x) and C im (t,x) are the concentrations in the mobile and immobile domains, respectively, in which x and t are the travel distance (here: column length) and time. v * is the pore velocity with is Darcy's velocity; i and K * are the hydraulic gradient and conductivity, respectively).
 H D is the hydrodynamic dispersion coefficient (with  L as longitudinal dispersion length and  p D as pore diffusion coefficient). Finally,   is the empirically determined DDMT rate constant. In natural porous media with more or less pronounced preferential flow paths,   can be composed of primary components  d attributed to molecular diffusion (e.g., van Genuchten & Wierenga, 1976), and secondary components  a attributed to other exchange processes such as slow advection (e.g., Gorelick et al., 2005;G. Liu et al., 2007). Both may contribute to the effective DDMT rate at varying proportions, that is, . One possible interpretation of the function  is, for instance, provided in both Gorelick et al. (2005) and G. Liu et al. (2007). The asterisks (*) in Equations 4 and 5, and this text section, mark the parameters related to processes that are assumed to depend, at a measurable level, on the experimental temperature. In contrast, other possible temperature effects, as for instance on the mechanical dispersion due to small variations of the grain diameters (see Section 2.2), are assumed to be negligible in this study.

Temperature Dependencies of Pore Diffusion Coefficient and DDMT Rate
Inversely linked by the tortuosity factor  , pore diffusion is directly related to molecular diffusion, that is, Parkhurst & Appelo, 2013). Hence, Equation 6, as adapted from Fogler (2006) and Parkhurst and Appelo (2013) and based on Equations 1 and 2, can be used to relate the pore diffusion intensity at a given temperature level to reference conditions, that is, Hereby, both the experimental (T exp ) and reference absolute temperature (T ref ) are measurable system variables. The diffusive component of the first-order proxy used in the DDMT may be assumed to be directly proportional to molecular diffusion, that is,   d m D (e.g., Coats & Smith, 1964;van Genuchten & Wierenga, 1976). This, in turn, implies that  d, and eventually also the effective mass transfer rate constant   , becomes a temperature-dependent parameter, and that with subscript indices analogous to Equation 6. The exponential f is an empirical, non-negative constant that allows to indirectly consider not only the additional contribution by slow advection, but also by other unknown or uncertain contributors such as geometrical parameters of the dual-domain system itself or hydrodynamic molecule radii. These can be subject to non-uniform temperature behavior. f is used in the evaluation procedure to quantify the deviation of the obtained  i temperature dependencies to a diffusion-only setup (with f set to 1).

Experimental Setup: Overview and Investigated Temperature Range
Eight 1-D column experiments, referred to as "E1" to "E8," were performed in this study to systematically assess the influence of experimental temperature on the transport behavior of tracers (see Table 1). In these experiments either (1) the experimental temperatures, (2) the employed tracer substance, or (3) both were varied, while all other experimental conditions remained constant.
The scope of the study was to investigate a temperature range that embraces both typical laboratory and field-scale settings. Groundwater temperatures typically reflect both natural (e.g., mean air tem-BINDER ET AL.

10.1029/2020WR029474
4 of 20 perature, geothermal gradients, river-groundwater interaction; Händel et al., 2013) and anthropogenic influences (e.g., urban subsurface structures, geothermal applications; Lund & Freeston, 2001) and, hence, may cover a rather wide temperature range. Therefore, temperatures ranging from 3°C to 30°C were considered (see Table S1). Besides these bounding values, two additional temperature levels, that is, 10°C and 20°C, were incorporated. The experimental series was not performed in a monotonously increasing or decreasing temperature order, but in an intentionally mixed order. This allowed to distinguish between temperature-related modifications and any possible trends associated with progressing experimental time. Furthermore, the initial 10°C experiment was repeated at the mid and the end of the experimental series to check for consistency and reproducibility.
All transport experiments employed a stainless-steel column setup ( Figure 1a) of cylindrical shape, with an effective length of 100.6 cm between the inflow and outflow filter membranes and an inner diameter of 7.6 cm. The column length variations due to thermal expansion/ shrinking was found to be negligible, with a <0.05 cm difference for the investigated temperature range ( 1 steel 1.6 10 K being the thermal length expansion coefficient of stainless steel and ΔT = 27 K). The column was filled with a poorly graded sand under wet conditions and with repeated consolidation. Inner spiral-shaped blockades were used to minimize the risk of preferential flow alongside the column wall. The sand was an artificially mixed, homogeneous sediment provided by a regional distributor for construction materials (Co. SAXONIA Baustoffe). Wet-sieving and sedimentation analysis (Casagrande aerometer method) de-BINDER ET AL.  termined the following mass fractions: 0% gravel, ∼36% coarse sand, ∼50% medium sand, ∼8% fine sand, and ∼6% fines (silt, clay). The saturated hydraulic conductivity at 10°C is ∼8 × 10 −6 m s −1 (falling-head permeameter method). For the investigated temperature range between 3°C and 30°C, the relative volumetric expansion of silica grains between lowest and highest temperature was estimated to be less than SiO 2 2 3 1.7 10 K as thermal volume expansion coefficient of silica). This indicates that both the sediment's total porosity (n ≈ 31.9%, mean value calculated for the column based on bulk and material densities employing the gas pycnometer method) as well as its uniformity coefficient can be considered as constant parameters. The same column setup was used in all experiments, that is, the experiments were performed without repacking.
The column was placed into a high-precision thermostat cabinet (Figure 1a, Co. POL-EKO-APARATURA), which allowed for an accurate temperature control across the entire temperature range, with fluctuations not exceeding ± 0.1 K. To avoid temperature-related density effects, the cabinet was equipped with a fan-regulated air distribution system to ensure a uniform temperature for the entire column. A schematic overview is provided in Figure 1b.

Experimental Setup: Procedure
KBr was applied in the experiments E1, E3, and E5 to E8 (Table 1) due to the known conservative behavior of the monovalent bromide anion (e.g., Levy & Chambers, 1987;Mastrocicco et al., 2011) and due to its very low background concentration. The EC of the injected water was ∼3.0 mS/cm (at 25°C) and the tracer-to-background difference, that is, ΔEC was ∼2.7 mS/cm. In the experiments E2 and E4 (Table 1), 2 H 2 O, that is, heavy water, which is generally considered as an almost perfect tracer (e.g., Becker & Coplen, 2001;Koeniger et al., 2010), was applied with an isotopic offset,  2 Δ H, of ∼94‰ between ambient (−59.6‰ vs. VSMOW) and injected water (+34.6‰ vs. VSMOW). When converted to absolute deuterium isotope concentrations (Becker & Coplen, 2001), this corresponds to a difference of ∼3.26 mg L −1 (tracer: 36.02 mg L −1 , background: 32.76 mg L −1 ). Both tracer solutions were prepared as single batches for the entire experimental series. Tracer stability of both KBr and 2 H 2 O was confirmed in complimentary batch experiments that were conducted at three temperature levels (10°C, 20°C, and 30°C) and with a 48-h sediment contact (Hahnewald, 2018).
For the propulsion of the tracers, a high-precision/low-pulsation peristaltic pump equipped with long-life PVC tubes (both Co. ISMATEC) was employed to implement a flux-based, that is, 2nd type boundary condition (BC) at the inflow. The water flux, that is, the pumping rate was set to a constant value of 0.382 g/ min, which was periodically (at least once per day) checked by determining the change of the water mass in the collection container. This pumping rate resulted in consistent Darcy flow rates of ∼12.1 cm/d for all experiments. Although small fluctuations were observed, differences in flow velocities typically remained within a very small range and the uncertainty was estimated to be significantly less than 1%. During temperature re-adjustments (e.g., from 10°C to 30°C), the pump was iteratively set to the new experimental condition until a stable flow was achieved. The outflow was designed as a 1st type, that is, constant head boundary condition. The employed 2nd-/1st-type combination ensured that temperature-related hydraulic conductivity changes only affected the hydraulic gradient between in and outflow, but not q Darcy directly. An alternative 1st-/1st-type approach with pre-defined head differences, that is, without a pump would have resulted in temperature-dependent variations of q Darcy (see Text S3 for an example).
Each tracer injection was performed for a precise period of 24 h, whereby only one tracer substance was injected per experiment to avoid detrimental chemical interactions that could have biased the tests otherwise. Furthermore, the tracer solutions were degassed to minimize outgassing and, hence, to avoid air entrapment within the sediment. After each tracer application, the column was continuously flushed with unlabeled and degassed tap water. The employed tap water was stored in a large, closed water barrel. It was stepwise transferred to gas-/water-dense, multilayered composite material pouches (Co. Tesseraux Spezialverpackungen). The tracers were stored in the same type of pouch, located outside of the thermostat cabinet due to limited space. However, the residence time within the inflow tube was sufficiently long to ensure a full temperature adjustment. This was checked for selected experiments by manual temperature measure-BINDER ET AL.
10.1029/2020WR029474 6 of 20 ments at the column inflow. One refill of the tap water barrel was necessary in-between the experiments E6 and E7, whereby the difference in the EC values of the two batches was less than 5 µS/cm. EC was used as proxy indicator for the salt tracer breakthrough. For this, a calibrated flow-through EC sensor, connected to a mobile data logger (both Co. WTW), was directly attached to the effluent end of the column. This allowed for a controlled monitoring at 5 min logging intervals. The EC signal, while being cheap and easy to measure, is a widely employed and standardized proxy indicator for saline tracers in environmental tracing and monitoring applications. The measurement procedure for EC is quickly adoptable to changing experimental conditions (e.g., Vienken et al., 2017).
In the experiments E1 to E4, discrete sampling was additionally performed at the column outflow end to determine either anion and cation concentrations as well as pH values, or deuterium isotope ratios (δ 2 H), respectively (Table 1). Sampling was carried out manually for ∼6.5 days, whereby each sampling round required ∼40 min. The interval between sampling rounds was adjusted to the observed breakthrough behavior with more frequent intervals during the peak (∼3-5 h pause between sampling events) and less frequent sampling intervals during the initial phase and during the tailing phase (pause of ∼11 -13 h, i.e., twice per day). Ion concentrations in the samples were measured using ion chromatography (Co. DIONEX) for anions and atom adsorption spectroscopy (Co. Varian) for cations, while δ 2 H was determined by using a high-precision isotope ratio mass spectrometer (IRMS) system (MAT 253, Co. Thermo Fisher Scientific GmbH) coupled to a high-temperature pyrolysis furnace (Co. HEKAtech GmbH). Temperature-dependent offsets in the EC measurements, among others due to changes of the employed sensor's cell constant, were directly addressed by standardized temperature compensation as described in the German guideline DIN EN 27888 (DIN, 1993). This procedure included a non-linear conversion from the measured EC value to a reference temperature (25°C) value to take the temperature-dependent relationship between ionic concentration and EC into consideration (DIN, 1993).
Note that the injection of the degassed tap water was not stopped after the end of the breakthrough monitoring, but continued until the start of the next experiment. Here, the average pause between experiments employing the same tracer substance was ∼2.5 total pore volumes (PV). This was done to ensure an ionic re-equilibration and, therefore, to provide consistent initial conditions for all experiments.

Data Evaluation Concept
For a straightforward comparison of the experimentally obtained BTCs and for the assessment of the effect of temperature change on the effective mass transfer rate constant, two dimensionless parameters were defined following Toride et al. (1993): the beta factor (Equation 8) describing the pore distribution and a first-order Damkoehler number (Equation 9), which characterizes the relation between the advective and DDTM timescales.
In general, the importance of considering DDMT-related effects typically increases with rising Da values, that is, with increasing mass transfer rates and with decreasing  values as the size and, thus, the storage and reaction capacity of the immobile domain increases. In this study, x was considered as a constant value, representing the length of the column.
Numerical solute transport modeling was used to analyze the experimentally obtained BTCs and to evaluate the hypothesized temperature dependence of . Our evaluation was performed for the effective DDMT rate constant as a whole, not specifically the diffusive part alone. The adjustable transport parameters ,  L, and  (respectively Da) were inversely determined using (1) (Harbaugh, 2005) and MT3DMS for flow and solute transport simulations, respectively (Tables S2 and S3). Subsequently, PEST (Doherty & Hunt, 2010) was employed to determine optimized parameter values for the adjustable transport parameters by minimizing the sum of the squared residuals between simulated and corresponding observed values (e.g., concentrations).
Four different calibration and, thus, data evaluation approaches (A1 to A4) were considered (Table S4). First (A1), the BTCs were separately evaluated employing the advection dispersion equation (ADE). Calibration parameters were  (by varying n m in Equation 8) and  L, with both being separately determined for each temperature level, while all effective mass transfer rates  i for temperature levels i were set to zero by definition. Next, the single-rate DDMT model was detachedly applied to each BTC data set (A2; additional calibration parameter group is , with separately determined  i), thereby assuming that the experiments were entirely independent from each other. In contrast, in A3 and A4 the BTCs were jointly evaluated to represent the hypothesis that both  and  L are temperature-independent and, hence, shared transport parameters valid for all temperature levels. Here,  i values were still considered to be separate, temperature-dependent calibration parameters. For A3, the ratios between the  i variations were not pre-defined, that is, the empirical constant f in Equation 7 was adjustable. In A4 the ratios between the  i variations were pre-defined by setting the empirical f value to "1" while assuming a temperature-dependence analogous to the known temperature dependence of diffusion coefficients. This final step A4 was implemented to evaluate through a separate model run if pre-defining changes in  i (assuming a diffusion-dominated mass transfer) could be used as additional constraint for solving the inverse problem. For all approaches A1 to A4, both the total porosity n (averaged value is 31.9% as derived in Section 2.2) and the pore diffusion coefficients D p,i were considered as constants with pre-defined values to reduce unnecessary ambiguity affecting the inverse modeling procedure. Hereby, D m and, hence, also D p values (with   p m D D n, i.e.,  approximated by n −1 ) were defined based on the water's ionic composition during the breakthrough peak considering the respective temperature level (Equation 6). It should be emphasized that self-diffusion (water molecules in water) was assumed for the isotope 2 H instead of tracer diffusion (as done for EC and bromide). Note that eight separate PEST runs with significantly varying starting conditions were performed for each approach in order to manually mimic a global optimization strategy. Figure 2 shows the EC-based BTCs that were obtained in the 24-h KBr application experiments, that is, E1, E3, as well as E5 to E8. Non-normalized concentration-time scatterplots are shown in Figures 2a and 2b, and a semi-logarithmic probability net graph is shown in Figure 2c. The probability net graph, in which any tracer transport unaffected by non-equilibrium processes will result in a quasi-straight line, allows for a better illustration of the first tracer arrival. It also allows for a preliminary parameter estimation for the basic transport parameters ( and  L), based on the graphical method proposed by Nitsche (1981;Text S2). Finally, a double-log view with power law (C ∼ t −m with m denoting the slope; see, e.g., Pedretti et al., 2013, for details) is provided in Figure 2d, which allows for an in-detail investigation of the late-time behavior following the peak. In Figure 2d, additional BTC data is presented after 2.6 total PVs (for 30°C until ∼6 total PVs).

Continuous EC-Based Observation Data
In all cases, the first breakthrough at the column's effluent end was observed after ∼1 day, that is, ∼0.4 exchanged total PVs, while the peak was recorded shortly after ∼0.8 total PVs (∼2 days). The shape of the BTCs points toward a small to medium contribution of the immobile region to the overall transport behavior, given that the post-peak tailing extends up to ∼3 total PVs (∼1 week). It is clearly evident that a change in the experimental temperature has a systematic effect on the shape of the BTCs. In particular, the EC data of the four different temperature levels result in a fanned BTC array (Figure 2a). It is apparent that higher experimental temperatures are associated with decreased peak values. The observed peak values for EC range between 1,623 µS/cm (30°C, E3) and 1,775 µS/cm (3°C, E7). When normalized to the maximum signal difference of the tracer (∼2.7 mS/cm), the EC difference of around ∼150 µS/cm corresponds to a sig-BINDER ET AL.
The first arrival and peak times of the EC signal at the column outflow are almost the same for all investigated temperature levels. The observed deviations in the first arrival times (a few minutes) are still within the error margin of the defined versus the actual flow rate. In terms of the late-time behavior, higher/lower experimental temperatures lead to a shift of the cumulative curves to the right/left in the probability net graph (Figure 2c), respectively. The results of the initial parameter analysis, employing the method after Nitsche (1981), changed as well (based on regression lines of the equilibrium ADE in Figure 2c). Between the temperature levels, estimated values range from ∼60% to ∼75% for  and from ∼0.5 to ∼1.5 cm for  L.
BINDER ET AL.
10.1029/2020WR029474 9 of 20 The tailing analysis employing the double-log points toward a non-constant slope m in the power law equation. In particular, slope m is ∼3 shortly after the BTC peak, and subsequently decreases to values between 1.8 and 2.0 (at ∼2 total PVs, i.e., after ∼5 days). During the very late experimental time (between four and six total PVs, i.e., more than 2 weeks), m noticeably rises again, eventually approaching a value of 5.
All described observations provide a strong indication that the modifications to the BTC shape are related to the temperature-dependent changes in the diffusion coefficients (pore diffusion and DDMT). Although less pronounced compared to the differences in the peak values, there is a post-peak reversal of the BTC order (Figures 2a and 2b) during the phase when the tracer is slowly released from the immobile regions, further underpinning this hypothesis. Note also that the BTCs intercept after the peak at approximately the same experimental time (after ∼1.1 total PVs). The same phenomenon was previously observed by van Genuchten and Wierenga (1976), especially for  below ∼0.1 d −1 , while performing a model-based and temperature-independent parameter study for the DDMT approach.
The data sets obtained from the 10°C experiments were almost congruent with only small differences in the observed EC values, mainly in the form of a small offset caused by the refill procedure of the tap water reservoir (compare experiments E1 and E5 vs. E8). Hence, there was no temporal trend. Furthermore, the comparably small variations for the EC base level at different temperature levels -with a maximum ΔEC of less than 20 µS/cm -indicate that analytical uncertainties related to, among others, the employed temperature compensation are significantly smaller than the observed differences in the breakthrough behavior.

δ 2 H and Bromide
For the 10°C and the 30°C temperature levels, discretely sampled BTCs for Bras well as δ 2 H are displayed in Figures 3a and 3b, together with their analytical errors. In general, the BTC shapes qualitatively mimic the EC observations discussed above. In particular, the first arrival and peak times were very similar, and the same finding more or less also applies to the normalized peak values, which were in the same range as for EC. In contrast to the clearly visible BTC differences obtained from the EC data, the discrete sampling, however, resulted in less pronounced deviations between the two analyzed temperature levels. The hypothesized impact of the temperature variations is, nevertheless, still visible. This is in particular true during the rising branch of the BTC. Here, the observed concentrations and water fractions, respectively, of the 30°C BTCs were continuously below the values of the 10°C BTCs. For instance, the isotopic deviation is ∼5‰ at 0.8 PVs (a few hours before the peak); this equals around five times the long-term IRMS precision.
Compared to the EC data, the temporal density of the discrete observation data is much lower, which indirectly leads to a larger contribution by single measurement errors. Besides that, all δ 2 H and Br − data represent values time-averaged over the length of the sampling event, which took significantly longer than the near-real-time EC measurement. This visual effect is further analyzed by down-sampling and volume-averaging the EC data to mimic the same sampling schedule as employed for the discrete sampling. As shown in comparative scatter plots (Figures 3c and 3d), the normalized data pairs (EC vs. Br − , EC vs. δ 2 H) mainly follow a 1:1 line, except for some outliers which can be attributed to tracer-specific behavior. When converted to relative fractions (not shown as BTC), very small yet still measurable differences in the breakthrough peak value and in the intensity of the tailing can be observed. Therefore, we hypothesize that the smaller, 2 H-labeled water molecules are slightly more influenced by diffusion than the effectively larger bromide anions. This also indicates that the observed temperature dependence of the EC-based BTCs (as shown in Section 3.1) is most likely a combined effect, as the EC represents various species with their own species-dependent diffusion coefficients. This would, in turn, also affect the investigated exchange between mobile and immobile regions. All in all, the quality of match between the employed tracer signals is within the expected limits previously observed in lab-scale (e.g., Koeniger et al., 2010) and field-scale (e.g., Luhmann et al., 2012) experiments employing combinations of stable isotopes and other tracers such as dyes and salts.
BINDER ET AL.

Further Sampling Results
The analysis of the remaining major ions indicates that ion exchange, which mainly caused a release of calcium, magnesium, as well as sodium cations, occurred during the tracer breakthrough at both investigated temperature levels, while potassium was temporally bound to the sediment's surface, eventually causing a significant retardation of this cation. Aside from the phenomena already discussed for EC, 2 H, and Br − , only small temperature-related differences were observed for the anion and cation BTCs (Text S2, Figure S3). Also, note in this context that, among others, the ionic charge balance in the liquid phase needs to be maintained at any time of the migration process (electrostatic neutrality

General Comments on the Calibration
The PEST-estimated parameter values ( ,  L, and ) for evaluation approaches A1 to A4 as well as accompanying fitness measures R 2 and RMSE are listed in Tables 2 (EC-based BTCs) and 3 (bromide and deuterium measurements). In these plots, the estimated mass transfer rates are shown as dimensionless Damkoehler numbers    Da to allow for a direct assessment of the relative impact. Further information on the calibration procedure is given in the supporting information, including detailed comparisons between simulation results and observations ( Figures S4-S16). Table S6 lists standard deviations for all calibrated parameters as part of the manual global optimization strategy employed by PEST, and Table S7 shows dimensional values for the mass transfer coefficient. In general, the optimization results highlight the importance of considering temperature effects on the estimated transport parameters. Observed particularities are explained in Sections 4.1.2 (EC) and 4.1.3 (bromide and deuterium).

Analysis of the Calibration Results -EC
Employing the ADE model (A1) leads to a comparably wide calibration parameter range with a difference of  Δ ≈ 7.6% and  Δ L ≈ 2.1 cm between 30°C and 3°C experiments. Qualitatively, the results follow a trend that can be anticipated from the BTC shapes (Figures 2a and 2c). However, the estimated  values are slightly higher and the  L values are much higher than the ones obtained from the graphical method ( Figure 2c). Furthermore, as the ADE does not account for non-equilibrium processes (e.g., Coats & Smith, 1964;Hunt et al., 2014), the residuals remain comparably large (Table 2; Figure S4). By employing BINDER ET AL. Note. The values given in parentheses (except for the 10°C value) denote the ratios between the respective averaged value and the reference at 10°C. a transport parameters were jointly assigned to all four BTCs/temperatures and fitted together. b value calculated based on pre-defined ratio to 10°C value. c initial tracer-to-background difference, that is, the reference value is ∼2.7 mS/cm. the DDMT model and, hence, by adding  (respectively Da) as a third adjustable parameter (A2), the residuals between simulation results and experimental data decrease significantly. Furthermore, the obtained values of  now differ only slightly, with  Δ ≈ 2.4%. Recalculated to porosities, n m and n im lie in the range between 21.3% and 22.1% and between 9.8% and 10.6%, respectively. Almost the same applies for the longitudinal dispersivity: on average, the estimated  L values are three times smaller than the corresponding values in A1, while  Δ L ≈ 0.9 cm also showed smaller differences. This is in accordance with previous findings (e.g., Coats & Smith, 1964). By employing approach A2, both  and  L decrease systematically with rising temperatures, especially the latter. This is the opposite behavior of that seen for the A1 results, which is also not immediately plausible given that lower peak values combined with wider BTC shapes are typically associated with increasing dispersivity values. In this scenario, the DDMT coefficient varies significantly, with inversely determined values for  i ranging between 6.5 × 10 −7 s −1 and 11.6 × 10 −7 s −1 (Da ≈ 0.47 to 0.83), correlating with experimental temperatures. For the weakly coupled parameter estimation approach (A3), the estimated  and  L represent approximately the average of the values obtained by the independent BTC simulation (A2), while the bandwidth of  i is slightly reduced (7.2 × 10 −7 s −1 to 10.5 × 10 −7 s −1 , i.e., Da ≈ 0.51 to 0.75). However, the change of this value still follows the same systematic behavior. Similar to A2, the computed residuals are again very small. Finally, if the empirical constant in Equation 7 is set to f = 1 (strongly coupled simulation A4), both  and  L exceed the range given by A3 and their estimated values are slightly higher than the estimates derived from the 3°C experiment (A2). Here, values estimated for the DDMT approaches range between 4.9 × 10 −7 s −1 and 10.7 × 10 −7 s −1 (Da ≈ 0.35 to 0.77), with slightly BINDER ET AL.
4.0 × 10 −1 ± 4.9 × 10 −2 4.9 × 10 −1 ± 6.3 × 10 −2 3.3 × 10 −1 ± 4.0 × 10 −2 30°C -6.9 × 10 −1 ± 1.0 × 10 −1 (172%) 6.0 × 10 −1 ± 6.8 × 10 −2 (122%) Note. The values given in parentheses behind the respective 30°C value denote the ratios to the references at 10°C. a transport parameters were jointly assigned to both temperatures and fitted together, but separately for bromide and deuterium data sets. b transport parameters were jointly assigned to both temperatures and fitted together, but separately for bromide and deuterium data sets. c value calculated based on pre-defined ratio to 10°C value. d initial tracer-to-background differences, that is, the reference values are ∼1,630 mg/L (Br) and ∼3.25 mg/L ( 2 H). e absolute isotope concentration is used here instead of the Delta notation. increased residuals. However, the residuals still remain below the ADE-based results (approach A1) after optimization. The decreased fitting quality (and also the increased residuals as shown in the supporting information) may indicate that for the investigated columns the temperature dependency of the overall mass transfer process cannot be solely attributed to diffusion.

Analysis of the Calibration Results -Bromide and Deuterium
The calibration results for bromide and deuterium (  Arrhenius plots for EC, bromide and deuterium showing  values calibrated with approaches A2, A3, and A4 (i.e., theoretical relation with f = 1). D p is included to function as reference having the same slope (with equidistant axis for  and D p ).
tions ranging from ∼0% to ∼7% for  (mean: ∼3%), from ∼0 to ∼4 cm for  L (mean: ∼1.2 cm) and from ∼0.1 to ∼0.3 for Da (mean: ∼0.2). Hereby, the deviations obtained for the A1 approach were more pronounced; and the deviations obtained for deuterium increased compared to those for bromide. Furthermore, the less dense temporal density of samples is clearly associated with a higher parameter uncertainty, especially with respect to Da.

Theoretical versus Observed Values
The inversely determined  i values showing variations between the different temperature levels were less pronounced than expected when compared to variations anticipated on the base of an analogous temperature dependence, as it is known for pure molecular diffusion. This observed discrepancy is illustrated in Figure 4a and further analyzed in Figures 4b-4d showing Arrhenius plots in which the natural logarithm of a measured parameter (here: D p and ) is plotted against 1/T. The y-axis intercept corresponds to the ln of the pre-exponential factor A E in the Arrhenius equation (see Stumm & Morgan, 1981), while the line's slope equals −E A /(N A × k B ). Here, E A is the activation energy and N A ≈ 6.02214076 × 10 23 mol −1 is Avogadro's constant.
In the approaches A2 and A3, the ratios between the  i (e.g.,    (Figures 4b-4d, with parameters in Table 4) for the investigated temperature range, but consistently shows a reduced slope compared to D p and  with pre-defined ratios (A4).
Overall, this is an indicator for a combined effect of diffusive and non-diffusive processes. Additional mass transfer attributed to slow advection within the porous media (G. Liu et al., 2007), which in turn may be less affected by the temperature change than its diffusive counterpart, is a possible explanation. Other explanations are more related to hydrochemical aspects. For instance, as the measured EC integrates the behavior of both conservative and reactive constituents (e.g., Mastrocicco et al., 2011;Singha et al., 2011), reactive processes may contribute in a non-linear fashion to the temperature-induced EC BTC modifications. Furthermore, variations of the electrochemical potential induced by the temperature change may affect the diffusive behavior of both the cations and anions (e.g., Rolle et al., 2018). Aside from biases caused by these hydrochemical processes, the observed attenuation of the temperature dependence may also be related to the physical characteristics of the porous medium itself such as, for instance, changing pore-scale geometrical parameters.

Conclusions
Several previous studies (e.g., C. Liu et al., 2008;Maloszewski et al., 1999) performed tracer experiments with multiple tracers in order to elucidate the impact of their varying diffusion coefficients on observed breakthrough characteristics. Following a similar approach, the present study investigated the transport of potassium bromide and deuterium oxide in a column filled with poorly graded sand and a small fraction of fines. Advancing from the previously reported studies, our experiments also investigated the influence of temperature variations. Together with the model-based interpretation, these experiments were specifically aimed at the identification and quantification of the temperature dependence of the mass transfer processes and how this would affect the solute breakthrough behavior in mildly heterogeneous porous media. Largely confirming our original hypothesis (Equation 7), both the experimental observations and the numerical modeling results suggest a significant influence of the ambient temperature on the transport of tracers. In all optimization approaches employing DDMT-based models, the effective mass transfer between mobile and immobile domains seems to follow an Arrhenius relationship, that is, the mass transfer rate constant increases with increasing experimental temperatures. However, our analysis points toward a contribution of additional processes that are not directly related to pure molecular BINDER ET AL.

10.1029/2020WR029474
15 of 20 diffusion, while showing a less pronounced dependency on temperature, and thus weakening the overall temperature effect.
While laboratory-obtained DDMT rates themselves can generally not be directly used for field-scale applications due to ubiquitous (up-) scaling effects, transferring reaction parameters appears to be more appropriate, as discussed, for example, in Greskowiak et al. (2010), C. Liu et al. (2008), and Li et al. (2006). In such cases, column studies and comparable laboratory-scale experiments employed for the determination of these parameters should ideally be performed under climate-regulated conditions at temperatures close to those present at the studied field sites to allow for a stringent and consistent chain of evidence. Additionally, this implies that experimental temperatures of complementary experiments (e.g., a reactive column experiment accompanied by a batch reactor study to determine decay coefficients) should not vary. For experimental setups comparable to the present study (1-D, similar ratios between advective and diffusive fluxes), we suggest keeping potential temperature variations within the same experimental series below a limit of ∼5 K. This "rule-of-thumb" recommendation is based on the observed parameter uncertainties of ∼10% ( ,  L) and ∼20% (), respectively, as determined from the EC data inversion employing the temperature-independent DDMT (approach A2). In case of larger variations, the effect of temperature, especially on the mass transfer, should be included in the BTC evaluation as demonstrated here. In this context, we also recommended to always report experimentally obtained mass transfer coefficients and reactive parameters together with the respective experimental temperatures. This is particularly true when EC is employed as a proxy for breakthrough detection, as additional biases introduced by, among others, temperature compensation algorithms must be considered.
Future investigations should expand the current results by additional experiments with a broader selection of sediment types and tracers as well as a wider range of flow velocities to further generalize this work, as the number of experimental data sets provided by the present study is still somewhat limited. These future experiments should be accompanied by systematic parameter studies that include reactive transport modeling to explore which processes compete with the DDMT process in terms of temperature dependency. Finally, the optimization methods used in this study (especially A3) may be exploited in laboratory-scale experiments to obtain DDMT coefficients with reduced uncertainties, as the combined observation data of BINDER ET AL. Note. Values were determined by regression analysis using the Arrhenius plot.

Table 4 Arrhenius Activation Energies E A and Accompanying Pre-Exponential Factors A E
multiple BTCs obtained at different temperature levels can then be used to constrain the inverse modeling procedure more tightly.

Data Availability Statement
Further information related to temperature bandwidths, model usage, and model results as well as to chemical analysis is provided in the Supporting Information S1. Digital versions of the data sets (primary experimental data, best-fit simulation results), shown in Figures 2a, 3a, 3b, and S2-S15, were uploaded to a permanent data repository (Figshare, see https://doi. org/10.6084/m9.figshare.11316791). Fichtner, Thomas Krause, Alireza Arab, and Sharif Ibne Ibrahim for many helpful scientific discussions. The authors would also like to thank the various members of the Department Environmental Informatics (UFZ, Leipzig) for useful suggestions and the staff members of the Institute of Groundwater Management (TU Dresden, Dresden) for their enduring technical support. Peter Kohl and René Hädicke are thanked for their contributions to preliminary experiments. Finally, we would also thank all of the anonymous reviewers who have provided valuable and constructive feedback. This study was funded by the German Research Foundation (grants 275519227, 402833446, and 185284906) and by the German Federal Environmental Foundation (grant 20015/363 personal to co-author P. Stock). Travel funding for a short-term research stay at the University of Western Australia (great!ipid4all program, grant 2016_55) was provided by the Graduate Academy of TU Dresden and by the German Academic Exchange Service (DAAD). The authors declare no conflict of interest. Open access funding enabled and organized by Projekt DEAL.