## 2 Materials and Methods
### Derjaguin-Landau-Vervey-Overbeek (DLVO) Interaction Energy Calculation  
To confirm the details of _n_ZnO transport behaviors in the presence/absence of SRHA, through the column experiments, the total interaction energy of NPs-NPs and quartz-NPs according to the distance was calculated using the classic DLVO theory. As shown in Equation (1), the DLVO interaction energy was calculated as the sum of retarded van der Waals energies (\(\Phi_{VDW}\)) and electrostatic double layer energies (\(\Phi_{EDL}\)). \[\Phi_{Total}=\Phi_{VDW}+\Phi_{EDL} \tag{1}\]  
In \(\Phi_{VDW}\), a sphere-sphere equation was used for NPs-NPs, and a sphere-plate equation for NPs-sand. Each equation is expressed as Equations (2) and (3) [40]. \[\Phi_{VDW-SS}=-\frac{A_{131}a_{p1}a_{p2}}{6h(a_{p1}+a_{p2})}\left[1-\frac{5.32 h}{\lambda}\ln\left(1+\frac{\lambda}{5.32h}\right)\right]^{-1} \tag{2}\]  
\[\Phi_{VDW-SP}=-\left\{\frac{A_{131}a_{p}}{6h}\right\}\left[1+\left(\frac{14h} {\lambda}\right)\right]^{-1} \tag{3}\]  
where \(a_{p1}\) and \(a_{p2}\) in Equation (2) refer to the radii of the two interacting spherical NPs, whereas \(a_{p1}\) in Equation (3) refers to the radius of the NPs, h represents the distance between the two NPs (Equation (2)) or the distance between NPs and sand surface (Equation (3)),and \(A_{131}\) is the Hamaker constant. In the NPs-water-NPs system, 1.23 \(\times\) 10\({}^{-20}\) J was used in the absence of SRHA, and 6.09 \(\times\) 10\({}^{-22}\) J in the presence of SRHA. \(A_{132}\) is the Hamaker constant in the NPs-water-collector system; it was calculated at 6.96 \(\times\) 10\({}^{-21}\) J in the absence of SRHA, and 1.54 \(\times\) 10\({}^{-21}\) J in the presence of SRHA. The combined Hamaker constant NPs-water-NPs system (\(A_{131}\)) and NPs-water-sand system (\(A_{132}\)) were calculated through Equations (4) and (5), respectively. \[A_{132}=(\sqrt{A_{11}-A_{33}})(\sqrt{A_{22}-A_{33}}) \tag{4}\]  
\[A_{131}=(\sqrt{A_{11}-A_{33}})^{2} \tag{5}\]  
\(A_{11}\) is the Hamaker constant for \(n\)ZnO, and 9.21 \(\times\) 10\({}^{-20}\) J was used [41]. In the presence of SRHA, 4.71 \(\times\) 10\({}^{-20}\) J was used for the Hamaker constant of \(n\)ZnO with adsorbed SRHA [42]. \(A_{22}\) was the Hamaker constant for silica, and 6.50 \(\times\) 10\({}^{-20}\) J was used, from [43]. \(A_{33}\) was the Hamaker constant for water, and 3.70 \(\times\) 10\({}^{-20}\) J was used same as \(A_{22}\), considering [43]. \(\lambda\) was the characteristic wavelength adsorbed in the interaction, and 100 nm was typically used. Initial particle sizes measured prior to transport experiments were used in Equations (2) and (3). Adsorption of SRHA on the surface of nanoparticles could stabilize the nanoparticles in suspension and alter particles size of nanoparticles, obtaining a different \(\alpha_{\rm p}\) used in Equation (3). To calculate \(\Phi_{EDL}\), NPs-NPs was assumed to have had a sphere-sphere form, and NPs-sand, a sphere-plate form. Each equation is shown as Equations (6) and (7) [44]. \[\Phi_{EDL-\rm{SS}}=\frac{2\pi a_{p1}a_{p2}n_{\infty}kT}{(a_{p1}+a_{p2})\kappa^{ 2}}(\phi^{2}{}_{p1}+\phi^{2}{}_{p2})\Bigg{\{}\frac{2\phi_{p1}\phi_{p2}}{\phi^{ 2}{}_{p2}+\phi^{2}{}_{p2}}\ln\left[\frac{1+\exp(-xh)}{1-\exp(-xh)}\right]+\ln [1-\exp(-2\kappa h)]\Bigg{\}} \tag{6}\]  
\[\Phi_{EDL-SP}=\pi\epsilon_{0}\epsilon_{r}a_{p}\bigg{\{}2\psi_{p}\psi_{c}\ln \left[\frac{1+\exp(-\kappa h)}{1-\exp(-xh)}\right]+(\psi^{2}{}_{p}+\psi^{2}{}_ {c})\ln[1-\exp(-2\kappa h)]\bigg{\}} \tag{7}\]  
\[\kappa=\sqrt{\frac{e^{2}\sum n_{i0}z_{i}^{2}}{\epsilon_{0}\epsilon_{r}kT}} \tag{8}\]  
\(n_{\infty}\) denotes the bulk number density of ions (N m\({}^{-3}\), where N is the number), \(k\) is the Boltzmann constant (J K\({}^{-1}\)), \(T\) the absolute temperature of the system (K), \(\phi_{p1}\) and \(\phi_{p2}\) are the zeta potentials of NPs 1 and 2 (\(\varphi=\mathrm{ze}\psi/kT\)), and \(\psi_{p}\) and \(\psi_{c}\) are the zeta potentials of the nanoparticle and sand (V), respectively.