Materials and Methods

Samples were collected from the Calhoun Critical Zone Observatory (CCZO) located in Union County, South Carolina, USA, as part of a multidisciplinary effort (supplementary Fig. S1). Related data sets for the CCZO are available and include listings for seasonal air and soil temperatures, climatic properties, LIDAR maps, groundwater and ground gas fluxes, photographs, soil properties, vegetative covers, and streamflow fluxes (criticalzone.org/calhoun/data/datasets/). In summary, the region’s annual precipitation averages 127 cm and the mean annual temperature is 15.7°C. The soils are composed of Cataula series (fine, kaolinitic, thermic oxyaquic Kanhapludults). Sites included in this study are pits located in experimental research watersheds that were excavated in 2016 using a backhoe to expose profiles down to ~3 m, which were sampled intensively. Deeper samples were extracted from the pits by hand auger to depths of ~8 m (Table 1). Cultivated plots in research watershed 1 (R1C) have been cultivated continuously with cover crops and occasionally amended with dolostone since 1930. Two pits excavated in R1C were considered in this study, with sites R1C2 and R1C3 located higher and lower, respectively, on the landscape. Research watershed 7 includes two pairs of forested plots (~15 m diameter) where all the trees have been identified, measured, and cored for age determinations. The forested plots were chosen with a hardwood plot that has been forested continuously, and a pine plot that was previously cultivated in each pair. Pits were dug in both sets of plots with the pit in pine plot #2 being the focus of this study (R7P2). Continuous sampling was conducted and bulk chemical and mineralogic analyses were performed for the profiles for all these sites, as well as four pits from two other watersheds; however, only a subset of shallow (<1 m) and deep (>6 m) samples was used from R1 and R7 as the focus for this current study (Table 1).

Table 1 Locations and properties of samples collected at the Calhoun CZO

More detailed site descriptions and data for the above and related samples are available at criticalzone.org/calhoun/data/datasets/

Samples were sieved to remove the >63 μm fraction and dispersed using a Branson Sonifier Cell Disruptor 350 (Branson Sonic Power Company, Danbury, Connecticut, USA) in a solution of 38 g of Na-hexametaphosphate (Alfa Aeser, Ward Hill, Massachusetts, USA) and 8 g of Na-carbonate (Baker Chemical Co., Phillipsburg, New Jersey, USA) per liter of deionized water. The clay fraction was separated to the <2 μm (equivalent spherical diameter) from the <63 μm fraction using centrifugation (Schroeder Reference Schroeder2018). All samples were considered to be Na-saturated after this treatment. Approximately 1 g each of K-saturated and Mg-saturated sample was prepared by exchanging with 1.0 M and 0.1 M KCl (Fisher Chemical, Fairlawn, New Jersey, USA) and MgCl₂ (Acros, Morris Plains, New Jersey, USA) solutions, respectively. Exchanges were repeated by centrifugation and solution renewal to ensure full saturation. After rinsing of excess salt with deionized water, clay slurries were sedimented and air dried on 11 cm² glass petrographic slides to ensure infinite thickness (0.1 g/cm²). Each slide was scanned over the range 2–32°2θ (0.01°2θ step and 0.1 s per step) using a Bruker D8 advance diffractometer (Bruker, Karlsruhe, Germany) with CoKα radiation (35 kV, 40 mA), goniometer radius of 21.7 cm, primary sollar slits, 0.6 mm scatter slit, Fe-CoKβ filter, 2.5° receiving slit, and a Lynx-Eye position sensitive detector. Data were collected for each sample in the states of air-dried (AD), ethylene glycol-saturated (EG) (24 h at 20°C), and heated overnight at 110, 350, and 550°C. This multi-specimen method was similar to that used effectively by Lanson et al. (Reference Lanson, Sakharov, Claret and Drits2009) on illite-smectites. Only the K- and Mg-saturated scans in the EG and 110°C state were considered for the comparison with calculated patterns, because the AD samples expressed multiple 0-water, 1-water, and 2-water hydration states (see discussion below).

Samples contained gibbsite (4.84 Å) and/or goethite (4.15 Å), which in each case allowed for correction of minor sample displacement errors (<±0.05 mm). Experimental data were Kα₂ stripped using Bruker Eva ^® software (Version 4.2.0.14) and exported in .xy format. A linear background was subtracted from each data set using the minimum count per second value and then imported as the experimental data into NEWMOD2 ^® software. NEWMOD2 was used to model the mixed-layer oriented XRD patterns (Reynolds Reference Reynolds, Brindley and Brown1980, Reference Reynolds1985; Yuan & Bish Reference Yuan and Bish2010). Tables S1a–f contain all NEWMOD2 model parameters used in the fitting for all XRD patterns.

A forward modeling approach was used, whereby parameters and abundances were adjusted to minimize the differences between experimental and model patterns. Parameters adjusted included: the d spacing of layer types, mean defect broadening distance, large number (N) of coherent scattering domain (low N = 3), ordering scheme (Reichweite), the percentage of layer types, and type abundance of exchangeable or fixed interlayer cation in 2:1 structures, dioctahedral or trioctahedral layer type, and abundance of octahedral iron in 2:1 structures. As discussed by Austin et al. (Reference Austin, Perry, Richter and Schroeder2018), using a dioctahedral structure is well suited to simulate the diffraction characteristics of structural octahedral Fe in biotite that has undergone oxidation from the ferrous to ferric state. The oxidation of Fe reduces symmetry of the trioctahedral structure and makes the octahedral site more dioctahedral-like in terms of XRD phenomena, (i.e. shortening of the b lattice parameter), which has been supported independently by far-infrared studies (Diaz et al. Reference Diaz, Robert, Schroeder and Prost2010). The layer charge must compensate for this loss by expelling positive interlayer cations (i.e. potassium). If the unit layer charge decreases from 1.0 to 0.75 then the structure becomes ‘illite-like.’ Hence, for the purpose of this study and modeling, the term ‘illite’ was used, following the convention of Barre et al. (Reference Barre, Velde and Abbadie2007a). This is not meant to be the same as authigenic illite formed during burial diagenesis or degraded muscovite (Schroeder Reference Schroeder1992).

The authors recognize that NEWMOD2 accommodates only two-layer types in model calculations and that, in nature, three-layer type mixed-layer clays or combinations of the same layer types in different hydration states may be present (Dumon & Van Ranst Reference Dumon and Van Ranst2016). Improvements to the current approach could be made by modeling of 0-water, 1-water, and 2-water layer types in 2:1 clay structures, particularly for XRD patterns in the AD state. The results, however, indicated that reasonably good simulations were possible in spite of this limitation and the sample treatment differences were robust enough to justify using the two-layer model. Using the fewest variables in a model is beneficial for the eventual application of clay mineral quantification in predictive capacities, such as work by Hillier & Butler (Reference Hillier and Butler2018) who predicted the extractable potassium properties of soils from XRD data. Given the reasonable goodness of fit results and the small visual difference seen between both experimental and model data, the two-layer model approach (i.e. NEWMOD2) was used.

Fits were assessed using Goodness of Fit (G), which was calculated using the method described by Toby (Reference Toby2006) using the equations:

(1) $R_{\mathrm{wp}}=\sqrt{\frac{{\sum }_i{w}_i{\left(y_{c,i}-y_{o,i}\right)}^2}{{\sum }_i{w}_i{\left(y_{o.i}\right)}^2}}$

(2) $R_{\exp }=\sqrt{\frac{N}{{\sum }_i{w}_i{\left(y_{o,i}\right)}^2}}$

(3) $G=\frac{R_{\mathrm{wp}}}{R_{\exp }}$

where y represents intensity values, w represents the weight (w = 1/σ²[y _o,i], σ²[y _o,i] = y _o,i), the subscript c indicates calculated counts, and the subscript o indicates measured counts of i°2θ (Austin et al. Reference Austin, Perry, Richter and Schroeder2018).

The shorthand syntax suggested by Schroeder (Reference Schroeder2018) was employed where the binary mixed-layer system is labeled ABXXRY, where A = smaller d-spacing layer type; B = larger d-spacing layer type; XX = percentage of layer A, and RY = Reichweite ordering scheme when Y = 0 R is random ordering, and Y = 1; R is nearest neighboring layer only dependence. For example, KS70R0 indicates kaolinite-smectite occurring with 70% kaolinite layer types, 30% smectite layer types, and random ordering. If a layer type is repeated, then the structure is considered a single structure (e.g. KK = pure kaolinite). In some cases, modeling the same sample exposed to different saturation and solvation states resulted in solutions where the layer types were the same, but the abundance of layer A vs B varied slightly. In these cases, the percentage of layer A (XX) was reported as a range of values. In other cases, the same sample exposed to different cation saturation and solvation states resulted in solutions where layer types and abundances changed. For NEWMOD2 modeling purposes, the small and large number of coherent scattering domains and the mean defect-free difference were varied in all cases (pure and mixed-layer phases) to minimize the difference between i _c and i _o.

To ensure accurate phase identification and quantification, XRD patterns of each sample under four different saturation and solvation states were modeled using NEWMOD2. During the modeling process, some theoretical considerations were required to ensure good fits that also made physical and chemical sense. Differences in the observed and modeled intensities at low angles (2–7°2θ; ~51–14 Å) could be only partially explained by the theoretical parameters accommodated in the NEWMOD2 model. The model mismatch for intensities at low angles sometimes resulted in modeled intensity being greater than the observed intensity. This was the case for all K-sat/EG-solvated samples and all Mg-sat/EG-solvated samples, except R12C3_600-650. In that case, a peak was observed at ~4°2θ (~25 Å) which was interpreted as illite-smectite. Three samples had observed intensities greater than modeled intensities (R1C2_13-80, Mg-sat/110°C & K-sat/110°C, and R1C3_600-650 K-sat/110°C).

This low-angle intensity of the XRD pattern is influenced in part by the Lp factor. Lp is calculated by

(4) $\mathrm{Lp}=\frac{1+{cos}^{2}2\mathrm{\theta \psi }}{2\theta }$

where ψ is the powder ring distribution factor (Reynolds Reference Reynolds1986). Briefly, this factor can have either a single crystal or random powder form. If the random powder form is used, the preferred orientation of the sample has a large effect on the resulting low-angle intensity. ψ is calculated using the standard deviation of the assumed Gaussian distribution of tilt angles of the random powder (σ*) and the size of the primary and secondary sollar slits (Reynolds Reference Reynolds1986). Variations in σ* cause the intensity at low angles to change, with higher larger σ* values (representing more random orientation) resulting in higher intensity.

NEWMOD2 allowed for parameter variations in the Lp function, which included σ* and both primary and secondary Soller slit sizes.. The Soller slit settings were kept constant throughout (6.6° primary and 2.5° receiving) for all calculations. To evaluate the differences in the observed intensities across treatments, the average intensities for all treatments of all samples from each complete profile (n = 34 for each treatment) were compared (Fig. 1). In both cation saturation cases, the EG-solvated samples had lower intensities in the low-angle region, thus indicating higher values for σ* and implying that these samples were more randomly oriented than the heated samples. As the samples were heated and the interlayer water and ethylene glycol were removed, the clay particles were assumed to behave in a manner that indicates more uniform orientation.

Fig. 1 Average XRD intensity (source = CoKα) of all samples grouped by treatment (n = 34 for each treatment). Color-coded regions represent one standard deviation with blended shades showing overlaps. Heated samples exhibited greater low-angle scattering compared to their EG-saturated states. K-saturated samples exhibited greater low-angle scattering compared to the Mg-saturated states. This response can be attributed partly to sample-sensitive changes related to orientation and scattering behavior (Reynolds Reference Reynolds1986). The increased low-angle scatter is consistent with a more powder-like Lorentz polarization (Lp) of the beam, whereas the decreased low-angle scatter, is consistent with a more single crystal-like Lp of the beam. Differences in low-angle scattering can also be attributed to mixed-layering, where the occurrence of quasi-super structures add intensity to the low-angle region

Reynolds (Reference Reynolds1986) determined that σ* = 12 is a good estimation for oriented clay slides prepared as described above. In an attempt to resolve the differences between modeled and observed intensities, σ* was adjusted during modeling. The calculation of G (Eq. 3) revealed negligible differences in quantitation as σ* was varied. Therefore, a value of 12 for σ* was used for consistency in the model fit for the range of 2–7°2θ (~51–14 Å) for samples in the Mg-sat/EG state (Table S1a–f).

In addition to differences observed in sample treatments, low-angle scattering was most pronounced in samples that included ordered mixed-layers and disordered smectitic mixed-layer minerals. This was attributed to mixed-layering, where the occurrence of quasi-super structures add intensity to the low-angle region. The lack of well defined higher order reflections that are expected in >R1 ordering make distinguishing between the effects of Lp and super structure low-angle scattering difficult. Although beyond the scope of this study, this may be resolved through a systematic study of particle orientations using rocking curve experiments (Reynolds Reference Reynolds1986).

The largest innovation to modeling mixed-layer clays presented in this paper is the use of NEWMOD2 for a systematic and quantitative comparison of mineral weight abundances (wt.%) across treatments. The method chosen for quantification required pre-assessment of various cation treatments and hydration states. This ensured that the suite of minerals purported to exist had a sound physical basis across all sample treatment patterns. The mixed-layer minerals and quantity of layer types were compared across the sample treatments using the following criteria to identify and assess discrepancies:

1. (1) The total number of kaolinite layers could not change with solvation, saturation, or heating.

2. (2) Expanding layer types could be collapsed by heating and K saturation, causing an apparent increase in illite-like layers.

3. (3) Changes in the d spacing of expandable layers could result in an apparent increase in kaolinite mixed-layers as long as no increase occurred in the abundance of kaolinite layers.

4. (4) Treatments were not expected to completely transform all 2:1 clay layers in a sample, e.g. some low-charge smectite layers may not completely collapse to 10 Å upon K saturation, resulting in an apparent increase in vermiculite layers.

5. (5) G values were minimized in compliance with the preceding rules.

In order to maintain compliance with the criterion (1), methods for the abundance calculations were limited to those model solutions that best balanced each layer type for all treatment/hydration conditions. This resulted in a dynamic tabulation of layer-type abundances (Fig. 2) using the NEWMOD2 solutions for each state. The numbered arrows are intended to give a sense of the quantity of layers that were transformed during the heating treatment. In some cases, the transformation showed more kaolinite layers than expected. For example, the K-sat/110°C treatment of R7P2_700-800 had fewer than expected kaolinite layers compared to the K-sat/EG treatment of R7P2_700-800. The XRD patterns for R7P2_700-800 K-sat treatments showed, upon heating, no increase in intensity or asymmetry on the low-angle side of the peak at 14°2θ (~7 Å), and a decrease in a broad area of higher intensity between 5 and 10°2θ (15–10 Å) centered around 7°2θ (~14 Å) (Fig. S7). This is interpreted as a collapse of expandable, mixed-layer kaolinite-vermiculite (~7 Å/~14 Å) to kaolinite-illite (~7 Å/~10 Å). This interpretation is consistent with the interpretation of the other treatments, though the best-fit model (both visually and with the lowest G) results in more kaolinite layers than expected. A similar response was observed in the R1C2_700-800 K-sat patterns, with a slight decrease in intensity on the low-angle side of the broad peak at 10°2θ (~10 Å) accompanied by an increase in both peak-width and intensity, on the low-angle side of the 14°2θ (~7 Å) peak (Fig. S3). Difficulties occurred in resolving NEWMOD2 model solutions that consistently provided mass balance to layer transformations for the K-sat samples.

Fig. 2 Quantitative representation of changes in modeled layer-type abundances based on changes in XRD patterns from different sample treatments (Table 2). Arrows depict the changes in layer types and abundances resulting from heat treatments. Values on the arrows quantitatively constrain the number of layer types transformed after heat treatment. Color codes show similar mixed-layer groups: kaolinite and kaolinite-smectite dominated (cool violet/blues), kaolinite-illite-like dominated (warm brown to yellows), and vermiculite dominated (greens). The shallow samples (left side) exhibited low abundances of expandable 2:1 layers in all profiles, resulting in relatively small but perceptible changes with treatments. Deep samples (right side) had more abundant/dynamic 2:1 illitic and expandable layers

In contrast, the Mg-sat/EG treatment for all samples was considered to be more consistent than all the K-sat sample treatments for mass balancing the layer types upon heating (Fig. 2). Differences in mixed-layer types manifested in K-sat states vs Mg-sat states were related to heterogeneous layer charges within the 2:1 layers and resultant layer heterogeneous dimensions (Lagaly Reference Lagaly1982; Gier et al. Reference Gier, Ottner and Johns1998). As noted by MacEwan & Wilson (Reference MacEwan and Wilson1980), structural contractions of K-sat states in 2:1 layers are variable in response to the amount of layer charge. The state of an octahedrally coordinated two-layer Mg-hydrate group [Mg(H₂O₆]²⁺ occupies a more uniform and consistent d spacing than K-sat states. Therefore, Mg-sat/EG was used for quantification.