1 Introduction

Concerns over excessive eutrophication have been increasingly reported (Fischer et al. 2017). As a consequence of high nutrient (especially phosphorus (P)) accumulation, eutrophication provokes nuisance algae blooms and “dead zoom” in aquatic ecosystem. Particular attentions have been made on appreciable subsurface P loss from agricultural land, which was well documented in coarse-textured soil, high weathered soil, and porous organic soil (Kleinman et al. 2015). Thus, it is an exigency to understand P transport in subsurface soil and subsequent leaching to surrounding water bodies.

Phosphorus transport is largely affected by soil structure and properties. Although metal oxides (e.g., magnetite, kyanite, and goethite) in soil enhance nutrients’ stability (Yaghi and Hartikainen 2013), macropores and preferential flow in subsurface soil provide direct pathways for P leaching (Bergström et al. 2015). An example is colloid-facilitated transport, a strong sorbent bonding P, would rapidly transfer P through the macropores to aquatic environment (Vendelboe et al. 2011).

To improve P acquisition, inorganic P fertilizers were reported to be co-applied with novel materials such as artificial humic substances (Yang et al. 2019), natural humic substances (Hua et al. 2008), and geochemical-modified sewage sludge (Zhang et al. 2020). Besides these materials, animal manure is always considered as an economical alternative to enhance P bioavailability (Shafqat and Pierzynski 2013). Manure is well known for its high proportion of nutrients and organic matter. But direct application of manure left severe environmental footprints, such as greenhouse gas emission and groundwater pollution (Chadwick et al. 2011). Thus, before application, composting is highly recommended to stabilize and sanitize the manure (Ren et al. 2018). The amount of humic acid (HA) in manure compost is a major indicator of product quality and biological maturity (Moral et al. 2009). And a high-quality manure compost exhibits an enrichment of soil HA and reduction of water-extracted organic acid in the soil biota (Bernal et al. 2009). It acts as a slow-release fertilizer because the humidification prevents nutrient loss from mineralization (Mulbry et al. 2005).

However, a significant role of HA on P loss has been largely underestimated, especially when mineral oxides extensively covered the soil. HA competes with nutrients for bonding receptors on mineral oxides through its functional groups (e.g., carboxylic acids, carbonyl, and phenolic hydroxyl groups) (Niu et al. 2011). In Fu et al.’s study, higher affinity was found between soil HA and synthetic goethite than P and goethite (Fu et al. 2013). Thus, soil retention capacity of nutrients would be deteriorated by soil HA and mineral oxide interaction. However, there is still a handful of studies that investigated P dynamic transport and retention under such HA and mineral-rich media.

In order to better understand the underlying mechanisms of low utilization efficiency and high leachability of P, we exanimated phosphate (PO43−) transport and adsorption in the soil environment by column and bench experiments. It is hypothesized that (1) HA (from soil natural material and decomposed manure) in the soil environment would simultaneously bind to inorganic soil surface; (2) the HA coating would eliminate the receptor sites on soil surface that were originally available for PO43−; (3) the application of HA could accelerate the P leaching in subsurface soil.

2 Materials and Methods

2.1 Sample Preparation

The silica sand (8 mesh) in this study was obtained from Fisher Scientific Inc. (Hampton, NH), which had a diameter of 2.36 mm, porosity of 0.4, and a surface area of 3.02 ± 0.02 m2 g−1 measured by BET gas adsorption (ASAP2010, Micromeritics, Norcross, GA). Silica sand was treated successively with sodium acetate, hydrogen peroxide, sodium dithionate, and sodium citrate to remove organic matter. Then it was saturated with Na+ using 1 M phosphate-buffered saline (pH 7.0). Treated silica sand had a bulk density of 1.48 g cm−3 and a cation exchange capacity of 15.2 cmol kg−1.

Goethite was prepared referring to previous study (Schwertmann et al. 1985). Briefly, after aging of 21 days at 25 °C, the suspension of the mixture of 1.0 M ferric nitrate and 1.0 M potassium hydroxide (1:9, v/v) was washed carefully with de-ionized water via centrifugation. The rinsed solid was then re-suspended in 0.4 M hydrochloride acid. After washed and dialyzed against de-ionized water, it was freeze-dried to obtain crystalline goethite. To coat silica sand, goethite was mixed with silica sand (1:10, 1:5, and 1:3, w/w) in 0.01 M sodium nitrate solution (pH 7.5) and shaken for 48 h. Coated silica sand was washed with 0.1 M sodium nitrate (pH 7.0) via centrifugation. After rinsed with de-ionized (DI) water, coated silica sand was oven-dried at 110 °C. In the study, X-ray photoelectron spectroscopy (XPS) analysis was used for determining the efficiency of goethite coating on silica sand. The X-ray tube setting was 60 kV and 55 mA. The quantification of goethite coating was determined by dissolving coated silica sand in nitric acid (95%) and hydrofluoric acid (40%) (2:1, v/v).

With the verbal permission from the farm owner, dairy manure was collected from a dairy farm located one mile south away from Mayo, Florida, for HA extraction (for research purposes only). The manure was mixed with DI water. And the pH of the mixed solution was adjusted to 12 with sodium hydroxide. After shaking for 1 h, the manure solution was filtered by a 0.45-μm filter. The pH of the filtrate was lowered with hydrochloride acid. Then the precipitates were collected and dried as extracted HA. HA was dissolved by de-ionized water to the concentrations of 0, 10, 20, and 40 mg L−1 for the following experiments. HA solutions have been used for both adsorption isotherm experiments and transport column experiments. Concentrations were quantified by UV-vis spectrometry (Agilent Technologies Cary 60) at a wavelength of 254 nm. The spectroscopic measurements were calibrated with a TOC analyzer (Shimadzu TOC 5000).

2.2 Adsorption Isotherm Experiments

Batch adsorption isotherms were conducted to examine P in the form of PO43− adsorption on goethite-coated silica sand following the procedures described in a previous study (Li et al. 2018). Estimating from previous studies, the commonly used range of HA is 0–20 mg L−1 (a typical soil bulk density = 1.6 g cm−3 was used for unit conversion) (Adani et al. 1998; Ferrara and Brunetti 2008; Tahir et al. 2011). In order to thoroughly estimate the effect of HA, we have also included the bad organic management using excessive HA (i.e., HA > 20 mg L−1). A series of 25 mL vials containing solutions (20 mL) with P (0–200 mg L−1) in the form of monopotassium phosphate (KH2PO4, pH = 7) and 4 g goethite-coated silica sand (HA was in range of 0–40 mg L−1) as well as blank controls (sealed with Teflon-lined screw caps) were agitated on a Wrist Action Shaker (Burrell Scientific, Model 75) to reach equilibrium. The suspension was then centrifuged at 12,000×g for 15 min, after which, PO43− concentration in the supernatant was measured by ionic chromatography. The amount of P adsorbed on goethite-coated silica sand was calculated based on the mass balance of the difference between the initial and equilibrated solutions. Three replicates were conducted for each sample.

In this research, P had Langmuir adsorption isotherms on goethite-coated silica sand. When equilibrium was reached, PO43− adsorption on goethite-coated silica sand was described as (Adachi et al. 1997):

$$ S=\frac{q_{\mathrm{max}}K{C}_{\mathrm{e}}}{1+K{C}_{\mathrm{e}}} $$
(1)

where S is the adsorbed P on goethite-coated silica sand; qmax is the maximum mass of P required to form a mono-layer coverage of the media surface; Ce is equilibrium P concentration; and K is the Langmuir adsorption constant.

2.3 Phosphorus Transport Column Experiments

Column experiments were conducted to study the impact of HA and mineral coating on the transport of anionic P in porous media. The column experiments were conducted under saturated conditions with the column (5 cm ID × 15 cm length) vertically oriented. The bottom of column was sealed with a custom fit to permit the flow of water and retain the silica sand. To pack the column with goethite-coated silica sand, the column was first filled with nanopure DI water to a height of 2 to 3 cm and the media were packed in the column through CO2 solvation to eliminate air pockets. A height of 2 to 3 cm of nanopure de-ionized water was maintained.

Before starting each experiment, approximately 10 pore volumes of nanopure de-ionized water (ionic strength adjusted with NaNO3 to 0.001 M) was eluted through the column to stabilize the column condition by a peristaltic pump (Masterflex, Cole-Parmer, Vernon Hills, IL). A conservative pulse tracer (chloride) breakthrough curve was generated separately before the introduction of P solution. For each experiment, 2 pore volumes of sample solution which contained K2HPO4 solution (50 mg L−1, pH = 7) and HA (0, 10, 20, and 40 mg L−1) were introduced into the column at a flow rate of 2.5 mL min−1. The column was then continuously flushed with nanopure DI water (ionic strength adjusted with NaNO3 to 0.001 M) until the background signal was detected by a fraction collector. HA and P in elution were measured by ionic chromatography and UV-vis spectrometry, separately. For each column experiment, it was duplicated for three times. And the inconsistency of breakthrough curves was within 5% (95% CI).

Phosphorus transport in goethite-coated silica sand was described by the equilibrium-kinetic two-site model (Toride et al. 1995):

$$ \upbeta R\frac{\partial {C}_1}{\partial T}=\frac{1}{\mathrm{P}}\frac{\partial^2{C}_1}{\partial {Z}^2}-\frac{\partial {C}_1}{\partial Z}-\upomega \left({C}_1-{C}_2\right)-{\upmu}_1{C}_1 $$
(2)
$$ \left(1-\upbeta \right)R\frac{\partial {C}_2}{\partial T}=\upomega \left({C}_1-{C}_2\right)-{\upmu}_2{C}_2 $$
(3)

where C1 and C2 are the dimensionless P concentration in the solution and on the media surface, respectively; β is the partition coefficient; R is the retardation factor (\( R=1+\frac{\uprho_{\mathrm{b}}{K}_{\mathrm{d}}}{\uptheta} \), where ρb is the media bulk density, Kd is the partition coefficient of P between the solution and media, and θ is the porosity); T is the dimensionless time; P is the Peclet number (\( \mathrm{P}=\frac{vL}{D} \), where v is the interstitial pore-water velocity, L is the length of column, and D is the dispersion coefficient); Z is the dimensionless axial coordinate; ω is the dimensionless mass transfer coefficient; and μ1 and μ2 are the first-order dimensionless deposition coefficient in the solution and on the media surface, respectively. For this research, it was assumed that the dominating P deposition occurred in the solution. Thus, μ2 was set to 0.

At the same time, a pulse-type boundary condition was used for the upper boundary and a zero gradient was assumed for the lower boundary (Yuan et al. 2017):

$$ {\left\lceil -\uptheta D\frac{\partial {C}_1}{\partial Z}+q{C}_1\right\rceil}_{Z=0}=\left\{\begin{array}{c}\mathrm{q}{C}_0,\kern0.5em \mathrm{at}\kern0.5em 0<\mathrm{t}\le {t}_0\\ {}0,\kern0.5em \mathrm{at}\kern0.5em \mathrm{t}>{t}_0\end{array}\right. $$
(4)
$$ {\left.\frac{\partial {C}_1}{\partial Z}\right|}_{Z=L}=0 $$
(5)

where q is the specific discharge (Darcian fluid flux); t0 is the duration of injection; and C0 is the initial PO43− concentration. Transport parameters in Eqs. (2) and (3) were obtained by fitting the experimentally obtained breakthrough data using an implicit, finite difference scheme of CXTFIT (Toride et al. 1995). All the parameters were optimized by minimizing the sum of squared differences between observed and fitted concentrations using the nonlinear least square method (Toride et al. 1995).

2.4 Statistical Analysis

The experimental data were processed and analyzed by one-way one-factor analysis of variance (ANOVA) using SPSS for Windows 14.0 (IBM, Armonk, NY). The corresponding Scheffé test was made to differentiate means with 95% confidence (p < 0.05) when the study of variance showed differences between means.

3 Results

3.1 Goethite-Coated Silica Sand Characterization

Different goethite/silica sand ratio (i.e., w/w ratio of 1:10, 1:5, and 1:3) led to variable goethite coating on silica sand. The corresponding goethite-coated silica sand samples are named as GS1, GS2, and GS3. The XPS results indicated an Fe peak at the energy level of 6.4 keV (Fig. 1). The peak was most pronounced for GS3 and least pronounced for GS1, with GS2 in between. Goethite coating of GS1, GS2, and GS3 was determined to be 7.5, 10.2, and 14.1 mg Fe g−1 sand, respectively. There was no obvious change in particle diameter and bulk density after goethite coating. However, after goethite coating, the surface area (S0) increased from 3.02 ± 0.02 to 4.91 ± 0.04, 6.78 ± 0.10, and 7.17 ± 0.84 m2 g−1, and the cation exchange capacity (CEC) decreased from 15.2 to 12.4, 9.75, and 7.21 cmol kg−1 for GS1, GS2, and GS3, respectively.

Fig. 1
figure 1

XPS analysis of pure silica sand, GS1 (goethite/silica sand w/w = 1/10), GS2 (goethite/silica sand w/w = 1/5), and GS3 (goethite/silica sand w/w = 1/3)

Full size image

3.2 Phosphorous Adsorption Isotherm on Goethite-Coated Silica Sand

Phosphorous adsorption isotherms on the goethite-coated silica sand in the absence and presence of HAs are demonstrated in Fig. 2. The adsorption of anionic P on GS1, GS2, and GS3 increased with the increase of the equilibrium concentration. P had the greatest adsorption on GS3, followed by GS2 and GS1, while HA limited the PO43− adsorption to different degrees. In this study, P exhibited Langmuir adsorption isotherms on goethite-coated silica sand. Great agreements between experimental data and simulated curves were indicated by the prediction coefficient R2 and significant coefficient p value noted in Fig. 2. Higher HA concentration generally led to a lower adsorption capacity of the media. This investigation agreed with previous assumption about the influence of HA.

Fig. 2
figure 2

Phosphorus adsorption isotherms data and Langmuir modeling curves on a GS1 (goethite/silica sand w/w = 1/10), b GS2 (goethite/silica sand w/w = 1/5), and c GS3 (goethite/silica sand w/w = 1/3) in the absence and presence of humic acid (HA), respectively (p < 0.01)

Full size image

3.3 Phosphorous Transport and Leachability in Goethite-Coated Silica Sand

Tracer (Cl−) transport was studied before P transport experiments in goethite-coated silica sand (i.e., G1, G2, and G3). Nearly all the input tracer was eluted from the column. Since tracer should not be retarded or retained in the media, the retardation factor was set to 1.0 (i.e., Kd = 0), and the deposition coefficients (i.e., μ1 = 0, μ2 = 0) were set to zero during the simulation. Dispersion coefficient D was thus determined and used for the following interpretation of P transport in the corresponding media.

Breakthrough curves of P transport in goethite-coated silica sand were characterized by a self-sharpening front, which became broader and diffuser at the elution limb (Fig. 3). The long-lasting tails indicated kinetic-controlled PO43− retention in the columns. Retardation was exhibited as the leading edge dominated by the reversible adsorption of PO43−. Decreased retardation was found when the concentration of HA was higher. In addition, negative charged P displayed irreversible adsorption or kinetic adsorption in goethite-coated silica sand, which was evidenced by the reduced mass recovery during the transport. P breakthrough curves were fitted well with the model. The accuracy of transport modeling was expressed by the sum of the squared differences which ranges from 0.92 to 0.99.

Fig. 3
figure 3

Phosphorus experimental and simulated transport breakthrough curves on a GS1(goethite/silica sand w/w = 1/10), b GS2 (goethite/silica sand w/w = 1/5), and c GS3 (goethite/silica sand w/w = 1/3) in the absence and presence of humic acid (HA), respectively (p < 0.01)

Full size image

4 Discussion

4.1 Influence of Goethite Coating on Phosphorus Retention

After coating with goethite, S0 and CEC changed linearly with the increased goethite dosage in the samples. CEC is an essential soil property in terms of nutrient management. Soil CEC is affected by soil properties such as parent materials, organic matters, and pH. Surface soil and wetland soil always exhibited large CEC because of intense reducing conditions and high dissolved organic matters (Guedes et al. 2015). In this study, reduced CEC in higher mineral-coated sand is caused by the neutralization of the charges through coating reaction (Scheidegger et al. 1993). Increased PO43− adsorbency and decreased CEC corresponded to the increasing dosage of Fe oxide, when HA = 0 mg L−1. At the same time, larger S0 was observed on GS2 and GS3. It has been well-estimated from previous study that the binding capacity of Fe oxides is proportional to its specific surface area (Gérard 2016). Thus, higher dosages of goethite in GS2 and GS3 contribute to the larger surface areas, which further provide abundant binding receptors for P adsorption to take place (Wang et al. 2012).

Parameter qmax, μ1, and R derived from Eqs. (1) and (2) indicate the goethite altered the soil properties (Fig. 4). Parameter qmax of GS3 appeared slightly ahead of GS2, but much higher than GS1 when HA = 0 mg L−1. The changes of qmax were correlated to the adsorption capacity changes of soil media. The phenomenon demonstrates that mineral coating efficiently trapped the mobile nutrients. μ1 demonstrates that P deposition on GS1 is less than GS2 and GS3. Parameter R derived from column experiments reveals a heavily retarded transport on GS3, which is also corresponding to its largest S0. Retardation increases the possibility of plant uptake during the P transport. When comparing simulated parameters (i.e., qmax, μ1, and R) of three samples, it is noticeable that higher goethite content significantly increased the retardation as well as the adsorption capacity. This conclusion is consistent with previous researches. Studies confirmed that the ability of soil constituents to bind PO43− is dominantly contributed to Fe/Al oxide and clay mineral (Gérard 2016; Zavaschi et al. 2020). In Bortoluzzi et al.’s study, Fe/Al sesquioxide distribution controls soil P sorption capacities regardless of the environmental gradient conditions (Bortoluzzi et al. 2015). Maluf et al. verified that a higher proportion of Fe2O3 and Al2O3 in Oxisol leads to a higher P adsorbency than the Entisol (Maluf et al. 2018).

Fig. 4
figure 4

a Retardation factor R as the function of humic acid, b deposition coefficient μ1 as the function of humic acid, and c maximum mass required to form a mono-layer coverage qmax as the function of humic acid, applying to GS1 (goethite/silica sand w/w = 1/10), GS2 (goethite/silica sand w/w = 1/5), and GS3 (goethite/silica sand w/w = 1/3), respectively (p < 0.01)

Full size image

4.2 Influence of HA on Phosphorus Transport and Loss in Soil

In the breakthrough curve (Fig. 3), retardation (resulted from the equilibrium adsorption) is exhibited by the delay of the curve, i.e., the ratio of solution velocity to that of PO43−. Retention and adsorption (resulted from kinetic adsorption) are presented by the reduced amplitudes, i.e., the ratio of the eluted PO43− to that of the input. HA causes a significant impairment on soil absorbency and retention, especially for GS3 (Fig. 3). The capacity of anionic P retention on GS3 was reduced ~ 70%. But the retardation was not severely infected. A clearer picture is revealed from Fig. 4 b and c; qmax and μ1 were all largely slumped after adding HA, while the changes of R were less pronounced after adding HA, as a flatter slope indicated in Fig. 4 a.

The possible interactions between HA, mineral oxides, and P have also been discussed in previous studies by mathematical models, filed investigations, and batch experiments (Fu et al. 2013; Gerke 2010; Wang et al. 2015). In this study, P affinity is weakened after HA-mineral binding formed, as the qmax was decreased by ~ 75% adding 40 mg L−1 of HA. Similarly, Wang et al. observed that less amount of PO43− was adsorbed on HA-coated ferrihydrite than the pure ferrihydrite (Wang et al. 2015). At the same time, less isoelectric point and less specific surface area were found on HA-coated ferrihydrite than pure ferrihydrite. The formation of bidentate ligands between PO43− and Fe oxide is the dominant mechanism estimated by attenuated total reflectance Fourier-transform infrared spectroscopy in their study.

The impacts of HA on P were various in previous studies, because of different soil properties, forms of nutrient, and even the extraction methods (Mumbach et al. 2020; Zavaschi et al. 2020). For example, Fu et al. verified higher P adsorption onto goethite was achieved with a lower pH either with the presence or absence of HA (Fu et al. 2013). Different absorbencies resulted from the formation of different ligands (i.e., monodentate and bidentate ligands) during pH variations. Borie et al. found that volcanic soil retained high humic substance as well as stable mineral-humus complexes (i.e., allophane, imogolite, Al- and Fe-humus complexes) (Borie et al. 2019). High organic content in volcanic soil reforms the accumulated P. More than 65% of soil P, which was previously claimed only bound to mineral, was linked to organic substances. da Rocha et al. concluded that P extracted from pasture was more recalcitrant to plant than that from agroforestry and forestry (Rocha Junior et al. 2018). Some studies determined that HA could be beneficial to the P bioavailability if managed properly. For example, HA significantly promoted the labile fractions of P from single superphosphate incubation in clayey soil and sandy soil (de Jesus Souza et al. 2019). HA-coated monoammonium phosphate (MAP) fertilizer retained more mobile P than pure MAP in mineral-rich soil (do Nascimento et al. 2018). This was because of the chelating function of HA. It combines cationic nutrients like Ca2+ and Fe3+, which simultaneously allows the inner P fertilizer to dissolve without precipitation or adsorption. Fink et al. derived the conclusion that organic matters only affect the P transport by retardation but not retention capacity of soil in the literature review (Fink et al. 2016). However, soil retention capacity and retardation were both changed after adding HA in this study (Fig. 3 and Fig. 4a). Similar conclusion was drawn, which indicated the mineral oxide coating has larger positive effect than the deterioration caused by HA on P adsorption (Yan et al. 2016). These various conclusions indicate that the largely unknown dynamic nutrient transport in the subsurface soil still existed.

5 Conclusions

The dynamic interaction among phosphorus (P), humic acid (HA), and mineral oxides has been clearly revealed in this study. Mineral oxides altered soil properties and then increased phosphate (PO43−) adsorption capacity, evidencing from the enlarged surface area and limited cation exchange capacity. However, the binding sites on mineral oxides preferred the ligand formation with HA to P. Minimal PO43− retention was observed in column experiments when HA was 40 mg L−1, indicating the majority of adsorption sites were covered by HA. The bioavailability of adsorbed P via monodentate and/or bidentate ligands is doubted, which reveals further investigation is needed. It can be concluded from this study that a HA-rich soil is not a guarantee for P-efficient utilization, especially when composted manure is applied. The economic benefits may not overcome the environmental footprint caused by such agricultural community. Therefore, it is an urgency for regulating HA-involved fertilizers.