Evidence of Hydrocarbon-Rich Fluid Interaction with Clays: Clay Mineralogy and Boron Isotope Data from Gulf of Cádiz Mud Volcano Sediments

Abstract

Clay dehydration at great depth generates fluids and overpressures in organic-rich sediments that can release isotopically light boron from mature organic matter, producing ¹⁰B-rich fluids. The B can be incorporated into the tetrahedral sites of authigenic illite during the illitization of smectite. Therefore, the crystal-chemical and geochemical characterization of illite, smectite or interlayered illite–smectite clay minerals can be an indicator of depth (temperature) and reactions with the basin fluids. The aim of this study was to determine the detailed clay mineralogy, B-content and isotopic composition in illite–smectite rich samples of mud volcanoes from the Gulf of Cádiz, in order to evaluate interactions of hydrocarbon-rich fluids with clays. Molecular modeling of the illite structure was performed, using electron density functional theory (DFT) methods to examine the phenomenon of B incorporation into illite at the atomic level. We found that it is energetically preferable for B to reside in the tetrahedral sites replacing Si atoms than in the interlayer of expandable clays. The B abundances in this study are high and consistent with previous results of B data on interstitial fluids, suggesting that hydrocarbon-related fluids approaching temperatures of methane generation (150 °C) are the likely source of B-rich illite in the studied samples.

1. Introduction

Mud volcanoes (MVs) are generated by extrusion activity involving the transport of clay-rich sediments, liquids and gases (mainly methane) from deeper regions to the surface [1,2,3,4,5,6,7,8]. In recent years, both the source of material and fluids have been the focus of research as they give us important information about the presence of hydrocarbon resources at depth or global methane fluxes to the atmosphere [9,10,11,12]. Extensive work has been done in the study of fluid sources and pathways in sedimentary basins, where a close relation exists between fluids and the nature of clays, as a result of clay dehydration at depth resulting in smectite illitization processes [13]. The illitization process generates fluids and overpressures at temperature ranges of ~80 to ~150 °C [14,15,16] and smectite is transformed to randomly interstratified (R0) illite–smectite minerals (I-S) and to more illitic ordered (R1–R3) I-S [14,15,16,17,18,19].

Boron is abundant in marine sediments [20] and sedimentary clay minerals illite/smectite (I-S) [21,22], which contain orders of magnitude more boron than other common diagenetic minerals (e.g., quartz, carbonates and feldspars). Boron is a highly mobile element, preferring aqueous phases to that of most minerals [23]. Thus, by understanding how the aqueous B is incorporated into clay minerals, important insights may be gained to the fluid and chemical dynamics of a sedimentary basin. To use this geochemical tool, one must be able to interpret the boron isotopic composition of paleofluids that were present in a basin at the time of clay mineral diagenesis. Boron is incorporated into the clay mineral structure in tetrahedral sites and can also be adsorbed to clay surfaces, including those in the clay interlayers [24,25]. B adsorption on clays causes a preferential ¹⁰B uptake in tetrahedral sites of the clay related to bond strength. There is a coordination change of B from trigonal in water (at pH < 7) to tetrahedral on the clay surfaces [26,27]. This fractionation of B isotopes between trigonal and tetrahedral coordination during fluid–rock interactions is temperature dependent and insensitive to mineral composition [24,28].

Several recent studies have highlighted the potential utility of B-isotope ratios as a tracer for fluid–rock interactions [25,26,29,30]. The adsorption of B on clay surfaces has been extensively studied [27,31,32], because it can be easily exchanged [26]. However, fixed-B is more useful for interpreting paleofluid B composition because the B-isotopic composition is fixed when B substitutes for Si as Si–O bonds are broken.

Thermal maturation of organic matter during the burial process produces oil, wet gas and dry gas (mainly methane). Numerous studies of light stable isotopes in clays (e.g., [32,33,34,35,36]) have shown that trace elements (N, B and Li) commonly found in I-S are associated with hydrocarbon-related fluids generated during the maturation of organic source rocks. In these studies, it was shown that the light isotope of each of these “heteroatoms” released from organic matter dominates the fluids, thus these trace elements are ideal tracers of organic inputs to pore fluids. Thus, organic matter can release considerable amounts of B, producing ¹⁰B-rich fluids [24]. Late-stage or deep diagenesis of clay minerals [37] coincides with the time/temperatures associated with organic maturation processes that lead to the expulsion and accumulation of hydrocarbons. Thus, ¹⁰B-rich fluids are a source of B that can be incorporated into the tetrahedral layers of illite during the process of illitization of smectite at depth. The authigenic illite preferentially incorporates ¹⁰B, thus the remaining fluids are relatively enriched in ¹¹B [24]. Therefore, the crystal-chemical and geochemical characterization of illite, smectite or interlayered illite–smectite (I-S) clay minerals can be an indicator of temperature and reactions with the basin fluids.

Molecular modeling is a useful tool for determining many aspects of minerals at atomistic scale helping the interpretation of many experimental phenomena related with minerals, especially clay minerals [38,39,40], including with borate anions [41]. In this work, Density Functional Theory (DFT) methods were used to obtain information about the incorporation of B into the clay mineral structure, for understanding the experimental results.

In the Gulf of Cádiz, extensive work has been done on the study of fluid sources and pathways [42,43,44,45,46], including basin-scale reactive-transport models [47]. These studies conclude that clay mineral dehydration during reaction of smectite to illite, from Mesozoic to Tertiary shale and marl units has been the major influence on fluid compositions in many of the Mud volcanoes [42]. Some samples from deep, hotter regions are associated with B and Li-rich fluid compositions, which have been associated with production of methane from organic rich sediments [48,49]. The chemical analysis of clay minerals and the study of the diagenetic evolution of these units is therefore of high interest in order to better characterize the fluid circulation system present in fluid venting areas such as the Gulf of Cádiz.

In this work, a detailed mineralogical characterization of samples coming from several mud volcanoes and the content and isotopic values of B in clay minerals were analyzed to determine the diagenetic evolution of clay minerals and the possible depth (temperature), origin and interaction of clay minerals present in the mud breccia with the methane bearing fluids.

2. Geological Setting

The Gulf of Cádiz is located at the front of the Gibraltar arc, the westernmost tectonic belt of the Alpine Mediterranean compressional system, which has formed in response to the convergence between the African and Eurasian plates. It has a complex geological history and has undergone several episodes of rifting, compression and strike-slip motion since the Triassic [50]. During Tortonian times (11.2–7.1 Ma), allochthonous units took place in the Gulf of Cádiz formed by the westward migration of the Alborán domain associated with the formation of the Betic–Rifian Arc. This Guadalquivir allochthonous unit [51] consists of a mixture of Triassic, Cretaceous, Paleogene and Neogene sedimentary units overlying the Paleozoic basement [50], and is responsible for diapirism of huge volumes of mud and salt of Triassic units and under-compacted Early–Middle Miocene plastic marls [50,51] (Figure 1).

Throughout this area, several methane gas-related seafloor structures have been identified, including mud volcanoes, areas of carbonate crusts and chimneys, gas pipes and mud mounds [51,52,53,54,55,56,57,58,59]. Gravity cores were collected during the ANASTASYA/01, MVSeis/08 and CHICA/11 cruises on board of Cornide de Saavedra and Hespérides. Samples come from short gravity cores (up 2 m long) taken on the top of 16 mud volcanoes of the Gulf of Cádiz.

Table 1. Geographical position of the MVs sampled in the framework of this study. All M3 samples correspond to a diapiric structure, not a mud volcano.

Mud Volcano	Samples	Water Depth (m)	Location
			Lat.	Long.
Diapir	M3 10-14	600	34°59′42.4″ N	6°49′50.4″ W
Mekness	M4 18-22	694	34°59′06.6″ N	7°04′21.6″ W
Almanzor	M8 58-62	1440	35°20′57.6″ N	7°30′40.4″ W
Pixie	M12 66-70	1639	35°20′13.2″ N	7°50′38.4″ W
MVSeis	M14 138-142	1611	35°23′42.6″ N	7°51′28.8″ W
Guadix	M25 34-38	1435	35°30′56.4″ N	7°32′45.0″ W
Cid	M26 50-54	1330	35°26′29.4″ N	7°29′04.2″ W
Boabdil	M27 18-22	1106	35°25′49.2″ N	7°10′45.0″ W
Gazul	M29 106-110	411	36°33′29.4″ N	6°56′06.0″ W
Albolote	C3 32-33	353	36°34′27.2″ N	6°52′46.6″ W
Tarsis	C4 36-37	550	36°29′17.4″ N	7°14′39.9″ W
Anastasya	C7 38-39, A2 48-50	457	36°31′20.9″ N	7°09′04.8″ W
Almazán	C9 84-85, A8 70-71	830	36°03′08.0″ N	7°20′01.6″ W
Aveiro	C10 27-28	1060	35°52′19.1″ N	7°26′15.2″ W
Faro	A14 26-30	795	36°05′31.8″ N	7°23′44.4″ W
Gades	A4 48-50	915	36°14′17.4″ N	7°37′01.2″ W

shows the location of the gravity cores of this study. Mud volcanoes are located in several zones in the Gulf of Cádiz (Figure 1) from 353 to 1639 m depth, more than 100 km apart from north to south. A detailed description of morphologies and processes can be found in [55,60]. Previous studies on the mineralogy of mud breccia [61] show that cores from the mud volcanoes are made of several units (hemipelagic and mud breccia) showing similar bulk mineralogical composition of quartz, clay minerals, scarce calcite and dolomite and pyrite, with slight differences in type of clay minerals between different mud volcanoes.

3. Experimental Methodology

3.1. X-ray Diffraction and Deconvolution

X-ray diffraction (XRD) patterns of oriented samples with a size fraction of <2 µm were obtained using a Bruker D8 Advance diffractometer, located at The University of Cádiz (Cádiz, Spain), with a graphite monochromator, operating at 40 kV and 40 mA using Cu-Kα radiation. Each sample was first washed with distilled water until the supernatant was chloride-free, sonicated and then the <2-µm fraction was separated by centrifugation [62]. Each suspension was smeared on glass slides and air dried in atmospheric conditions. The slides were then saturated with ethylene glycol at 80 °C for 24 h to ensure maximum saturation. XRD patterns were acquired on the oriented clay mounts in both air dried and ethylene glycol saturated state to determine the percent of illite in I-S [63]. To discriminate between detrital smectite and I-S mixed-layer phases, deconvolution of the patterns obtained from the oriented mount after glycolation were performed using the MacDiff 4.2.6 program (4.2.6, Johann Wolfgang Goethe-Universität, Frankfurt, Germany). The determination of the illite percentage (% illite) and type of order (Reichweite; R) in I-S was performed according to the position of XRD peaks 001/002 and 002/003 in the regions 8–11° 2θ and 14–19° 2θ, respectively [64,65]

3.2. Transmission Electron Microscopy

Grain morphology within the bulk and <2-µm fractions and quantitative chemical analyses by analytical electron microscopy (AEM) were obtained using a Philips CM20 transmission electron microscope (TEM) at the University of Granada (Granada, Spain). Powdered portions deposited on a holey C-coated Au grid were used to collect AEM spectra in scanning transmission electron microscopy (STEM) mode on areas of 200 Å× 1000 Å using a 70-Å diameter spot size. To check volatilization of light elements, analyses were taken at 15 and 40 s. The structural formulae of smectites, micas and interstratified I-S were calculated on the basis of 22 negative charges, i.e., O₁₀(OH)₂, adjusting the occupation of the octahedral sheet to 2 atoms per formula unit.

3.3. Secondary Ion Mass Spectrometry

Secondary ion mass spectrometry (SIMS) (Arizona State University SIMS Facility, Tempe, AZ, USA) was used to characterize the content and isotopic composition of B in the clay minerals. Analytical protocols for measurement of B contents and δ¹¹B values have been described elsewhere [49,66,67,68] and particularly for measurements of clay minerals by [24].

3.3.1. Sample Preparation

Boron is strongly adsorbed to the surfaces of clay minerals at room temperature, with a distribution coefficient >30 [68]. During burial, B can be found in two sites of clay minerals: exchangeable B in the interlayer and substituted B in the tetrahedral layer. Therefore, special preparations are required to separate exchangeable B from that held in the silicate framework before any isotope analysis of the structurally substituted B. First, the samples were treated with a 1N solution of mannitol a B-complexing polyhydric alcohol, which removes exterior surface B contamination [69], but not clay interlayer B. Samples were sonicated in an ultrasonic disaggregator, centrifuged at high speed to concentrate particles (or clusters of minerals), and then washed in triplicate in “B-free” deionized water filtered through Amberlite resin [33]. An aliquot of the mannitol-treated sample was mounted for isotope analysis of the total B (B _tetrahedral + B _interlayer) by drying a 5-mL suspension onto a one-inch (25-mm) diameter B-free glass slide. Several samples were placed on a single round B-free glass slide, including standards. Then, the measurement of total-B content was determined by SIMS using a calibration curve based on the counts of B (mass 11) relative to Si (mass 30). The calibration curve was measured on standard reference materials with known B-content [33].

The remaining clay was cation exchanged with 1 N NH₄Cl by standard procedures [62] to remove exchangeable B from the interlayer [70]. Samples were rinsed again in mannitol and then mounted for isotope analysis as above. These samples only contain B substituted in tetrahedral sites.

3.3.2. Boron Content and Isotope Analysis

A Cameca IMS 6f at Arizona State University (Tempe, AZ, USA) was used with a primary beam of mass-filtered ¹⁶O^– ions accelerated at 12.5 kV onto the sample held at 9 kV for a total impact energy of ~21.5 kV. Primary beam currents below 10 nA were used with beam diameters defocused to 40–60 µm. Positive secondary ions were accelerated away from the sample, and energy filtering (−75 V sample offset) was used for measurements of B-content [66]. No energy filtering was used for isotope ratio measurements.

B isotope ratios are reported in delta notation as:

δ¹¹B = [{(¹¹B/¹⁰Bsample)/(¹¹B/¹⁰Bstandard) − 1}*1000] − IMF

(1)

where the standard is NBS SRM 951, boric acid, with a ¹¹B/¹⁰B ratio of 4.0437 [71]. The instrumental mass fractionation (IMF) is determined by measuring a mineral standard on which the δ¹¹B is known. B-isotope analyses were calibrated by measuring clay mineral standard IMt-1 (Silver Hill Illite) from the Clay Minerals Repository (http://www.clays.org/sourceclays) that had been characterized by bulk thermal ionization mass spectrometry (TIMS) [24], with a δ¹¹B of −9 ± 0.6‰. The isotope ratio analyses averaged 50–100 cycles of measurements on each spot (depending on the B-content) and analytical errors were compared to predicted errors. Where analytical errors were >2 times predicted errors, the analysis was discarded. Multiple spots were analyzed on each sample and results were averaged. The internal standard was measured in between analyses of the unknowns to test for changes in IMF due to instrumental drift.

3.4. Computational Methodology

We created models of Al(OH)₃ and B(OH)₃ molecules enveloped in a hydrogen-bonding network of water molecules which simulates B- and Al-rich fluids that are present in the illitization process (at low pH). Two different models of illite structures were also created: one of them with tetrahedral Al (Al-illite) and the other one with B incorporated to the tetrahedral layer by replacing the Al (B-illite). The comparison of energies of these optimized components can show us whether the incorporation of B in the tetrahedral sites is energetically favorable. Furthermore, montmorillonite models were created to compare optimization energies between montmorillonite with B as hydroxide (B(OH)₄^–) in the interlayer and montmorillonite with B in the tetrahedral layer, replacing Si sites.

The electronic structure of Al(OH)₃ and B(OH)₃ molecules was studied by quantum chemical calculations with the Hartree–Fock approximation and the second-order Moeller–Plesset method for all electrons. A triple-ζ basis set with polarization functions was used for all atoms including H atoms (MP2/6-311G** level) as implemented in the Gaussian03 program package [72]. All geometries were fully optimized using the Berny analytical gradient method. No geometry constraint was applied to the molecules. Normal mode analyses were performed to the same level to confirm the nature of the various stationary points, finding only positive eigenvalues for minima.

Ab initio total energy calculations of the periodic illite crystal models and Al(OH)₃ and B(OH)₃ hydrated models were performed using density functional theory (DFT) methods implemented in the SIESTA program (version 3.0, Max Centre of Excellence, Modena, Italy) [73]. The generalized gradient approximation (GGA) was used with the Perdew–Burke–Ernzerhof (PBEsol) parameterization of the exchange-correlation function optimized for solids [74]. Core electrons were replaced by norm-conserving pseudopotentials [75]. Calculations were restricted to the Γ point in the irreducible wedge of the Brillouin zone. In all structures, the geometry of each atom was relaxed by means of conjugated gradient optimizations at constant experimental volume. In SIESTA, the basis sets are made of strictly localized numerical atomic orbitals (NAOs) with a localization cut-off radius corresponding to an energy shift of 270 meV. The basis sets used here are double-Z polarized (DZP) following the perturbative polarization scheme. This approach was successfully used in previous calculations on phyllosilicates [76] and hydrated systems [77,78,79]. A uniform mesh with appropriate plane-wave cut-off energy is used to represent the electron density, the local part of the pseudopotential, and the Hartree and exchange-correlation potentials. Total energy calculations were performed with cut-off energy values of 350 Ry. These conditions are consistent with previous studies with phyllosilicates [39,80].

Models

Based on previous works reporting quantum mechanical calculations [81,82], models of hydrated Al(OH)₃ and B(OH)₃ molecules were created, consisting of Al(OH)₃ or B(OH)₃ molecules optimized at MP2/6-311G** level, encaged in a cavity of a hydrogen-bonded network formed by 24 water molecules. Those models were also optimized using the DFT methodology implemented in the SIESTA program in the same conditions as the mineral structures.

Illite models were based on previous pyrophyllite models [38]. Pyrophyllite is a dioctahedral phyllosilicate [83] with a structure similar to illite, but without cation substitutions causing the layer charge on basal siloxane surfaces. The trans-vacant crystal form was used in all models [39]. To obtain a reasonable size illite model, a 4 × 2 × 1 supercell was generated. Two types of illite models were created, Al-illite and B-illite. The Al-illite model was generated from the supercell by replacing eight tetrahedrally coordinated Si⁴⁺ by Al³⁺, and four octahedral Al³⁺ were replaced by Mg²⁺. Layer charge is balanced by twelve K⁺ cations per supercell in the interlayer, resulting in a simulation cell composition of [K₁₂][Al₂₈Mg₄][Si₅₆Al₈]O₁₆₀(OH)₃₂. The B-illite model was created similar to Al-illite model, but replacing one of tetrahedral Al³⁺ by B³⁺, resulting in a simulation cell composition of [K₁₂] [Al₂₈Mg₄][Si₅₆A₇B₁]O₁₆₀(OH)₃₂.

In both cases, maximum dispersion of the substituted cations in the tetrahedral and octahedral sheets was made according to previous studies [40,84]. Initial lattice parameters of each 4 × 2 × 1 illite supercell are a = 21.14 Å, b = 18.35 Å, c = 9.79 Å; α = 91°, β = 100°, γ = 90°.

Montmorillonite models were created with a unit cell of a = 5.16 Å, b = 8.97 Å, c = 13.61 Å; α = 91.2°, β = 100.5°, γ = 89.6°, leaving enough space in the interlayer for avoiding additional variables related with interlayer complexes. Supercells of 2 × 2 × 1 were generated by replacing one octahedral Al³⁺ by Mg²⁺. Layer charge is balanced with one K⁺ cation per supercell. Two montmorillonite models were created: One with the salt K⁺ B(OH)₄^- in the interlayer, resulting in a simulation cell composition of [B(OH)₄^-][K₂][Al₁₅Mg₁][Si₃₂]O₈₀(OH)₁₆, and the other with B replacing one tetrahedral Si atom and Si(OH)₄ in the interlayer, resulting in a simulation cell composition of [Si(OH)₄][K₂] [Al₁₅Mg₁][Si₃₁B₁]O₈₀(OH)₁₆. In the illite and montmorillonite models the effect of the presence of water molecules can be considered similar in both cases with B and without B complex for the study on which is focused this work. Then, water molecules were not included to avoid additional computational effort and convergency problems found in our preliminary calculations.

4. Results

4.1. Clay Mineralogy

XRD and TEM analyses of samples were performed to characterize the clay minerals. Although preliminary XRD data indicated that samples were mainly smectite [61], the deconvolution (using MacDiff 4.2.6) of the pattern obtained from the oriented mount after glycolation in the regions 8–11° 2θ and 14–19° 2θ showed that detrital micas and mixed-layer I-S phases are present in addition to smectite (Figure 2). I-S from all samples have similar characteristics, presenting Reichweite (R) values corresponding to both the R0 (random ordering) and R1 (nearest neighbor) ordering stages (

Table 2. I-S peaks characteristic of the deconvoluted XRD patterns. Illite proportions were calculated following standard [64] procedures.

Sample	R	°2θ (001/002)/% Illite	°2θ (002/003)/% Illite	°Δ2θ/% Illite
M3 10-14	R0	10.31/10%	15.75/10%	5.44/10%
	R0	9.85/45%	16.25/50%	6.40/45%
	R1	9.17/80%	16.83/75%	7.66/75%
	Illite	8.78	17.61	8.93
M4 18-22	R0	10.29/10%	15.83/15%	5.54/15%
	R0	9.82/45%	16.16/45%	6.34/45%
	R1	9.18/80%	17.07/80%	7.89/80%
	Illite	8.83	17.51	8.68
M8 58-62	R0	10.32/10%	15.78/10%	5.46/10%
	R1	9.64/55%	16.61/60%	6.97/60%
	Illite	8.84	17.72	8.88
M12 66-70	R0	10.35/10%	15.82/10%	5.47/10%
	R0	9.99/40%	16.04/35%	6.05/35%
	R1	9.54/65%	16.70/60%	7.16/65%
	R1	9.30/75%	16.97/75%	7.67/75%
	Illite	8.89	17.73	8.84
M14 132-138	R0	10.22/15%	15.86/15%	5.64/15%
	R1	9.78/50%	16.48/55%	6.70/55%
	R1	9.36/75%	16.90/75%	7.54/75%
	Illite	8.82	17.56	8.74
M25 34-38	R1	9.72/55%	16.45/55%	6.73/55%
	R1	9.16/80%	17.04/80%	7.88/80%
	Illite	8.73	17.71	8.98
M26 50-54	R0	10.23/15%	15.82/15%	5.59/15%
	R1	9.74/55%	16.44/55%	6.70/55%
	Illite	8.79	17.65	8.86
M27 18-22	R0	10.17/25%	15.96/25%	5.79/25%
	R1	9.68/55%	16.40/55%	6.72/55%
	R1	9.22/75%	16.85/75%	7.63/75%
	Illite	8.83	17.53	8.70
M29 106-110	R0	10.05/30%	15.99/30%	5.94/30%
	R1	9.74/55%	16.32/55%	6.58/55%
	R1	9.19/80%	17.11/80%	7.92/80%
	Illite	8.86	17.82	8.96
C3 32-33	R0	10.30/10%	15.80/10%	5.50/10%
	R0	9.76/50%	16.28/50%	6.52/50%
	Illite	8.87	17.81	8.94
C4 36-37	R0	10.09/30%	15.97/25%	5.88/25%
	R1	9.81/45%	16.29/50%	6.48/50%
	Illite	8.84	17.74	8.9
C7 36-37	R0	10.06/30%	15.98/25%	5.92/30%
	R1	9.55/60%	16.58/60%	7.03/60%
	Illite	8.92	17.74	8.82
C8 23-24	R0	10.07/30%	15.99/30%	5.92/30%
	R1	9.68/55%	16.52/55%	6.84/55%
	Illite	8.97	17.69	8.72
C9 84-85	R0	10.28/15%	15.84/15%	5.56/15%
	R1	9.67/55%	16.46/55%	6.79/55%
	Illite	8.87	17.75	8.88
C10 27-28	R1	9.18/80%	17.06/80%	7.88/80%
	Illite	8.87	17.81	8.94
A2 48-50	R0	10.22/20%	15.83/20%	5.61/20%
	R0	9.81/50%	16.35/50%	6.54/50%
	R1	9.31/75%	16.92/75%	7.61/75%
	Illite	8.81	17.69	8.88
A4 48-50	R0	10.08/30%	15.98/25%	5.90/30%
	R1	9.66/55%	16.47/55%	6.81/55%
	R1	9.33/75%	17.04/80%	7.71/75%
	Illite	8.90	17.80	8.90
A8 110-111	R0	9.77/50%	16.19/45%	6.42/45%
	Illite	8.84	17.78	8.94
A14 26-30	R0	9.94/40%	16.12/40%	6.18/40%
	Illite	8.79	17.36	8.57

).

The deconvolution of the peaks corresponding to 001/002 and 002/003 reflections into small peaks allows estimating the proportion of illite in each I/S, because each d value corresponds to R-ordering type.

A detailed TEM study combining chemistry and the imagery of individual clay–mineral particles on the holey carbon grid showed the expected morphological differences between the smectite and I-S: smectite displayed typical aggregates and flakes of irregular and wavy outlines and I-S showed individual crystals with more euhedral shapes and clear outlines (Figure 3). These morphologies are consistent with previous TEM observations on smectite illitization [85]. Kaolin hexagonal particles and palygorskite fibers were also present in the samples in minor amounts (Figure 3).

Figure 4 shows compositional variation of individual clay crystals determined by analytical electron microscopy (AEM). Figure 4a is a plot of tetrahedral Al content vs. interlayer K, showing trends that relate to particle morphology as defined by [85]. Chemical ranges correlate with TEM observations, where smectite flakes, aggregates, round and polygonal particles are characterized. The substitution of Al by Fe in the octahedral sites provides a well-defined negative relationship between these two elements (Figure 4b).

A similar negative relationship is also observed between Al and Mg (Figure 4c). No significant trends were found between Si and Fe or Si and Mg (not shown). However, there is a poorly defined positive relationship between Fe and Mg (Figure 4d). There is a negative relationship between Al ^TOT and Si (Figure 4e), indicating the generation of tetrahedral charge. Figure 4f shows the relationship between layer charge and Si content, where mainly beidellitic samples have more than half of the layer charge produced by tetrahedral substitutions. However, montmorillonite samples can be also represented by larger amounts of Si.

The plot of data in MR3–2R3–3R2 diagram [86] (Figure 5) clearly shows different chemical ranges in the samples analyzed, from montmorillonite and beidellite to illite compositions. The presence of I-S phases is represented by multiple points with compositions between smectite and illite.

4.2. Boron Content and Isotopes

The bulk B-content of the clay minerals studied range from 97 to 146 ppm. Cation exchanged (xc) samples contain slightly less B than the bulk mannitol washed (mw) samples, due to the loss of interlayer B after cation exchange. This implies that most of the B in the clay is fixed in the tetrahedral layer (structurally bound) and only minor amounts of B were held in the I-S interlayer. Figure 6 shows that there are no spatial correlations between δ¹¹B values and mud volcanoes locations or bathymetric depth.

4.3. Computational Modeling

Geometry optimizations of Al-illite and B-illite models were performed at constant and variable volumes. The lattice parameters and the main geometrical features of the optimized structures are consistent with experimental values (

Table 3. Cell parameters and select interatomic distances of Al-illite and B-illite unit cells (distances in Å and angles in °).

Features	Exp a	Al-illite	B-Illite
a	5.22	5.24	5.24
b	9.02	9.10	9.08
c	10.07	10.09	10.08
d(001)	10.02	10.02	10.03
α	90.0	101.1	101.1
β	95.7	95.8	95.9
γ	90.0	89.9	89.9
d(Si-O)	1.65	1.67	1.67
d(B-O)			1.48
d(AlVI-O)	1.96	1.95	1.94
d(O-H)		0.97	0.97

^a XRD experimental data [91].

), with a mean basal d(001) value of 10.02 Å for Al-illite and 10.03 Å for B-illite. The main geometrical features of the hydrated Al(OH)₃ and B(OH)₃ models (Figure 6) are presented in

Table 4. Interatomic distances of hydrated Al(OH)₃ and B(OH)₃ models (distances in Å and angles in °).

Average Features
Hydrated Al(OH)3 Model		Hydrated B(OH)3 Model
d(Al-O)	1.75	d(B-O)	1.37
d(O-H)	1.01	d(O-H)	1.01
α(Al-O-H)	120.11	α(B-O-H)	110.83
d(Hw-Ow)	1.05	d(Hw-Ow)	1.05
d(Hw…Ow)	1.50	d(Hw…Ow)	1.55
d(HAl…Ow)	1.76	d(HB…Ow)	1.67

. This table shows average values being smaller bond distances in B(OH)₃ molecule than in Al(OH)₃. These bond lengths are consistent with B-O and Al-O distances observed in the B-illite model described in Table 4. In the B hydrated complexes, the structure is highly symmetric with the cation and O atoms in the same plane according to previous studies at gas phase [88], whereas the H atoms are twisted to a different plane and oriented towards the vicinal O atom of the same molecule. However, in the Al(OH)₃ hydrated complex, the Al cation forms a pyramidal form with the O atoms; similar non-planar configurations have been reported in Al hydrates previously [89]. There are strong H bonding interactions between the metal hydroxides and the surrounding water molecules. Hence, different B-O bond lengths are found d(B-O1) = 1.355 Å, d(B-O2) = 1.371 Å and d(B-O3) = 1.394 Å and consequently different BO-H bond lengths are found, d(O1-H) = 1.023 Å, d(O2-H) = 1.004 Å and d(O3-H) = 1.002 Å. This can be explained due to a strong H bond with one water molecule d(O1H…Ow) = 1.563 Å, d(O2H…Ow) = 1.633 Å, d(BO2…Hw) = 1.495 Å, d(O3H…Ow) = 1.795 Å, Some correlation can be observed: a stronger H bond, a longer O-H bond and a stronger B-O bond length is. A similar effect is observed in Al(OH)₃ with d(Al-O1) = 1.743 Å, d(Al-O2) = 1.787 Å and d(Al-O3) = 1.836 Å and also strong H bonds d(O1H…Ow) = 1.554 Å, d(O2H…Ow) = 1.640 Å and d(AlO2…Hw) = 1.262 Å, d(O3H…Ow) = 1.567 Å. During the optimizations, proton exchanges are observed by dissociation of water molecules, according to previous works on Al hydrates [90].

The optimization energies of the two illite models (Figure 7) are compared with the energies of the hydrated Al(OH)₃ and B(OH)₃ models by proposing the next reaction:

U_Al-illite + U_B(OH)3 ↔ _B-illite + U_Al(OH)3

(2)

where U is the internal energy of the system. Our results shows that UAl-illite+UB(OH)₃ is less stable (−105751.2580 eV) than UB-illite + UAl(OH)₃ (−105751.6327 eV) in 36.7 KJ/mol. This means that the fixation of B on the illite tetrahedral layer instead of Al atoms is energetically favorable.

The optimization energies of the two montmorillonite models (Figure 8) are also compared, showing that montmorillonite with tetrahedral B is more stable (−48414.1023 eV) than montmorillonite with interlayer B(OH)₄⁻ (−48413.1163 eV), again proving that the incorporation of B in the tetrahedral layer is a energetically favorable process. Furthermore, in the optimized model of montmorillonite with interlayer K⁺B(OH)₄⁻, it can be noticed that this salt is dissociated in B(OH)₃ and K⁺(OH)⁻. The B(OH)₃ forms a H bond with basal tetrahedral O atoms, d(O…H) = 1.75 Å, and with the OH anion d(O…H) = 1.68 Å. Our energetic result is consistent with previous empirical studies [92] that indicated B(OH)₃ dominates in the I-S interlayer space.

Further calculations will be performed exploring these phyllosilicate models with different moisture grade (water molecules) and several pressure conditions (sediments environments) to complete this study but they are out of the scope of the present work.

5. Discussion

5.1. Clay Mineralogy and Diagenetic Evolution of Deep Sediments in the Gulf of Cádiz

TEM images and XRD data show that smectite and interstratified illite–smectite (I-S) are present as well as other clay minerals such as detrital mica, kaolin plates and palygorskite. The chemical analyses conducted on individual particles, indicate that all smectite and I-S phases are dioctahedral, with Al as the main octahedral cation. Notably, Si from the smectite layers in the mixed-layered I-S decreases as illitization progresses (Figure 5). Focusing on the nature of interlayer cations, the compositional analysis shows a wide range of variations in the content of K, Na, Ca and Mg, but the dominant cation is K. There is a good correlation between the interlayer K content and tetrahedral Al (Figure 4a). From this plot, together with the observations of these morphologies by TEM (Figure 3), it can be inferred that the increase in the illitization produces morphological and chemical changes that increases K and tetrahedral Al in clay particles. This is supported by the percentage of illite in I-S and Reichweite (R) estimations (Table 2 and Figure 5).

Some smectite–illite transformation characteristics can be observed in the MR³⁺–2R³⁺–3R²⁺ ternary diagram (Figure 5) [86]. In this diagram, the chemical composition of clay particles is represented together with ideal compositions of the discrete phases involved in the transformation (montmorillonite, beidellite and illite). The dioctahedral character of the phases involved in the illitization process is represented by the distance of the data from the 3R²⁺ apex. Furthermore, for illite-rich phases, an increase in the layer charge and interlayer cations can be seen, represented by the trend of the data towards the MR³⁺ apex.

All the clay geochemical trends are related to prograde diagenetic changes, where the dominant clay minerals provide clues about the burial/thermal history of sedimentary basins [19]. The reduction of expandable layers or the smectite–illite transformation processes have been related to the evolution of petroleum systems, as illitization of smectite overlaps the oil window [93,94,95,96,97,98], as discussed by numerous papers addressing the potential for reaction through solid state transformation or dissolution/precipitation mechanisms [17,19,61,85]. Although the coexistence of different I-S phases in a sample can correlate with a prograde evolution in a diagenetic sequence, in this study, the physical mixture of clay minerals may come from different stratigraphic units as a result of fluid expulsion during emplacement of the mud volcano [13]. Nevertheless, morphological changes observed by TEM (polygonal to euhedral) point to similar mechanisms described previously [85] in a series of experimental hydrothermal conditions as seen by X-ray diffraction and TEM, where they suggest illitization mechanism driven by dissolution/crystallization processes.

MV clay mineral samples in this study, with δ¹¹B values ranging from +2.2 to 12.7 and an average of −2.2‰ fixed-B abundances, are relatively high, 82–145 ppm. Again, a precipitation process is the most probable mechanism to incorporate B in the tetrahedral layer of illite as described previously by the authors of [68,99], showing that, during diagenesis, as temperatures approach 120 °C, B-adsorption becomes negligible and substitution of Si by B occurs as illite forms [100]. Molecular models presented in this work are in agreement with this statement, as these calculations prove that it is energetically favorable for B to reside in the tetrahedral sites of illite. Although B-O bonds in the tetrahedral layer are shorter (1.48 Å) than Al-O and Si-O bonds (1.78 and 1.67 Å, respectively), the lattice parameters of B-illite models are similar to the Al-illite model and to the experimental illite values, meaning that the incorporation of B in tetrahedral sites has no effect on the crystal structure.

Boron geochemistry has been studied recently in different environments as an indicator of fluid circulation and diagenetic grade. A set of samples from different mud volcanoes around the world as indicators of progressive diagenesis show a good correlation between B contents and δ¹¹B isotopic values [101,102], although these studies did not carefully separate tetrahedral B from interlayer (trigonal) B. Besides B, the uptake of N in illite (as NH₄⁺ substituting for K⁺) increases during diagenesis as illitization proceeds and has been studied in hydrocarbon-bearing sedimentary basins suggesting a kerogen source for both elements [32,103,104,105]

The nature of clay and clasts present in the mud breccia of the Gulf of Cádiz mud volcanoes can be used to infer the possible depth of the underlying units [42,56]. In Yuma mud volcano in the Moroccan margin, more than 200 clasts from the mud breccia were studied [106], displaying a very complex mixture of material from the sedimentary successions. The reconstructed sedimentary succession showed sediments at least as old as Eocene, with the presence of several clayey units Miocene in age (Aquitanian and Tortonian).

The detailed clay mineral characterization and the chemical composition made in this study indicates that although I-S content can vary among the mud volcanoes, the clay mineralogy is similar to that found in Tertiary units (Miocene in age) common in Mediterranean Messinian sediments. The Messinian clay minerals taken in DSDP legs 13 and 42A contained large amounts of smectite and are constant throughout the Mediterranean Basin [107,108]. Similar results were found by [109] and [110] in the Lower Messinian in Sicily or by [111] in the Upper Cenomanian–Turonian sediments in the high Atlas in Morocco. In addition, clays (dioctahedral smectite and illite) are a common component in the Miocene–Pliocene lithostratigraphic formations (Gibraleon clays) of the lower Guadalquivir Basin [112,113,114].

Based on the δ¹¹B values of I-S in MVs of the Gulf of Cádiz (−12.7‰ to +5.3‰;

Table 5. Boron isotope analysis. mw, Mannitol washed samples; xc, NH₄Cl exchanged samples. Standard IMt-1 is analyzed to determine IMF for each analytical session. This value is subtracted from the delta value. SE is standard error of the average. PE is predicted error, which is the best possible error based on counting statistics.

Sample	Preparation	11/10	IMF	δ11B	SE ‰	PE ‰	n	B (ppm)
M3 10-14	mw	39.192	−27.8	−3	1.0	0.9	3	142
	xc	39.151	−27.8	−4	0.3	0.8	3	132
M4 18-22	mw	39.139	−27.8	−4.3	0.6	1.0	3	104
	xc	39.090	−27.8	−5.5	0.5	0.5	2	88
M8 58-62	mw	39.568	−27.8	6.3	2.0	1.4	3	96
	xc	39.361	−27.8	1.2	0.5	0.4	2	82
M12 66-70	mw	39.163	−27.8	−3.7	0.9	0.5	2	99
	xc	39.115	−27.8	−4.9	0.1	0.6	2	89
M14 138-142	mw	39.030	−27.8	−7	0.1	0.6	2	105
	xc	38.965	−27.8	−8.6	0.5	0.5	2	102
M25 34-38	mw	39.305	−24.6	−3.4	0.6	0.4	2	137
	xc	39.192	−24.6	−6.2	0.6	0.4	2	134
M26 50-54	mw	39.224	−24.6	−5.4	0.6	0.3	2	101
	xc	39.204	−24.6	−5.9	0.5	0.4	2	94
M27 18-22	mw	39.442	−24.6	0	0.7	0.3	2	146
	xc	39.236	−24.6	−5.1	0.8	0.4	2	141
M29 106-110	mw	39.337	−24.6	−2.6	0.5	0.3	2	132
	xc	39.171	−24.6	−6.7	0.6	0.3	2	122
C3 32-33	mw	39.685	−24.6	6	0.6	0.3	2	126
	xc	39.657	−24.6	5.3	0.8	0.5	2	115
C4 36-37	mw	39.644	−24.6	5	0.7	0.5	2	113
	xc	39.499	−24.6	1.4	0.6	0.4	2	103
C7 38-39	mw	39.309	−24.6	−3.3	0.6	0.4	2	120
	xc	39.058	−24.6	−9.5	0.6	0.4	2	113
C8 23-24	mw	39.596	−24.6	3.8	0.9	0.5	2	163
	xc	39.297	−24.6	−3.6	0.7	0.4	2	145
C9 84-85	mw	39.220	−31.8	1.7	0.5	0.4	2	97
	xc	39.050	−31.8	−2.5	0.5	0.4	2	96
C10 27-28	mw	39.001	−31.8	−3.7	0.6	0.4	2	137
	xc	38.892	−31.8	−6.4	0.7	0.4	2	121
A2 48-50	mw	38.908	−31.8	−6	0.7	0.4	2	139
	xc	38.840	−31.8	−7.7	0.8	0.4	2	138
A4 48-50	mw	38.913	−31.8	−5.9	0.8	0.4	2	129
	xc	38.638	−31.8	−12.7	0.8	0.4	2	117
A8 110-111	mw	39.446	−31.8	7.3	0.8	0.4	2	111
	xc	39.200	−31.8	1.2	0.7	0.4	2	110
A14 26-30	mw	39.628	−31.8	11.8	0.7	0.5	2	125
	xc	39.240	−31.8	2.2	0.7	0.4	2	122

), the average clays equilibrated in fluid with δ¹¹B < 10‰, (calculated using the mineral-water fractionation factor 1000 ln α _min-water = 3.28 − 10.35 (1000/T); [33], which is significantly more ¹⁰B-enriched than seawater (+39‰). We interpret the B-isotope composition of the I-S in MVs of the Gulf of Cádiz to result from illitization of the smectite-rich sediments, probably from Messinian sources over a range of temperature during organic maduration of primary source rocks (~80–150 °C). Hydrocarbon-related fluids generated at temperatures of methane production (~150 °C) are enriched in ¹⁰B [24,115], thus hydrocarbon related fluids are the likely source of isotopically light B in I-S in all the studied samples, from north to south and from the shelf to 1639 m depth. There are no mineralogical sources in these sediments that could be a source of such high concentrations of ¹⁰B-enriched fluids.

5.2. B Isotopes and Origin of Fluids in the Gulf of Cádiz

The affinity of the light isotope ¹⁰B for tetrahedral coordination and the heavy isotope ¹¹B for trigonal coordination was shown by [27]. Hence, during illitization, clay minerals will concentrate ¹⁰B in the process of crystallization under hydrothermal conditions [24]. It was shown [116] that kerogen in the Gulf of Mexico sedimentary basin oil source rocks have a B isotopic composition of −4‰ to +10‰. B-isotopic values from Gulf of Cadiz mud volcano samples are similar to those previously reported by [116] (Figure 9). However, mud volcano samples show higher dispersion of δ¹¹B data ranging from −7‰ to +11.8‰ for the bulk (mw) samples and from −12.7‰ to +5.3‰ for the cation exchanged (xc) clay samples (Table 5), perhaps reflecting mineralogical variations arising from mixed fluid sources, as expected in the processes of expulsion of fluids in a mud volcano. Table 5 shows that there are small differences in isotopic compositions between mannitol and cation exchanged samples indicating that the interlayer B is isotopically heavier than the tetrahedral-B.

Figure 10 represents B concentration and B isotope composition of MV from different locations. It can be seen that the pore fluids from MV clay minerals are ¹¹B-enriched compared to the clay minerals. While this may in part reflect the mineral–water fractionation of B, it also reflects the mixing of basinal water with seawater (+39.5‰). Using the B-isotope fractionation equation for pH < 7 [33], the fractionation between illite and water at 150 °C (methane generation temperature) is −1.2‰. Therefore, an illite precipitated with δ¹¹B −12.7‰ (Table 5, Sample A4 48-50 (xc)) would be in equilibrium with +8.4‰ water. This isotopically light δ¹¹B water composition, compared to seawater, is consistent with ¹⁰B enrichment from organic sources. Furthermore, Gulf of Cádiz MV samples show δ¹¹B compositions near values associated with source rock kerogen samples in other hydrocarbon-rich sedimentary basins (e.g., [99]). The higher δ11B values measured in some samples (up to +5.3‰) may represent more mixing with seawater or samples less influenced by hydrocarbon-related waters. Seawater contains only about 5 ppm B, thus the high B content of the mud volcanoes requires a high fluid:rock ratio, as expected for MVs. Furthermore, where the water pH was higher than ~8, B(OH)₄⁻ dominates the fluid, and, because ¹⁰B prefers tetrahedral coordination, there is little fractionation between B(OH)₄⁻ dominated fluid and illite [28].

These data are consistent with previous results shown by [43,45], based on Li and Sr isotope data from interstitial water, suggesting that the fluid geochemistry in MVs from the Gulf of Cádiz is influenced by overpressuring caused by clay dehydration at several kilometers depth. Temperatures approached 150 °C where thermogenic methane is produced and illitization ends. These clay and marl units underneath mud volcanoes in the Gulf of Cádiz, are Mesozoic and Tertiary in age (mainly Tortonian-Messinian) and have been considered as the oil source rock in the area [44]. A complex scenario of hydrocarbon fluid generation at depth and migration to the upper units, dehydration of clays and mixing with methane-rich fluids, and a later mixing with shallow gas was proposed by [121] to explain the distinct composition of fluids of this area. The detailed clay mineral characterization together with the B-isotopic composition is used in this study to provide additional information on the fluid geochemistry obtained in this area by the authors of [43,44,45,46,121,122]. These new data allow us to corroborate that illitization is an important process in the generation of fluids in MVs and can be one of the driving forces of mud volcanism in the area. Boron geochemistry is also relevant in elucidating both mineral processes and fluid origins, as previously proposed by the authors of [25,43,102]. This also is in agreement with other works that show the importance of the relationship of clays and petroleum in oil formation, migration, accumulation and storage [122,123].

6. Conclusions

Mud volcanoes in the marine environment usually involve a mixture of clays dragged up by fluids from underlying units of various depths. The detailed characterization of the clays, as well as B isotopic compositions contributes to our understanding of the geological model of the area. Those data show that mixed-layer illite–smectite phases and other clay minerals present in the mud volcano samples were derived from depths where temperatures were great enough to generate B from organic source rock. The illitization process occurs at temperatures close to oil generation [37]. During the illitization process, B released from kerogen, enriched in ¹⁰B, was incorporated into the tetrahedral layer of diagenetic illite. The Gulf of Cadiz MVs are dominated by minerals with high B-content and low δ¹¹B, suggesting that they formed at depth, in equilibrium with hydrocarbon-related fluids at temperatures hot enough to have generated methane that is associated with these MVs. This interpretation is supported by theoretical atomistic calculations demonstrating the preferred incorporation of ¹⁰B in the tetrahedral sheet rather than in the interlayer space of the I-S. From an oil industry point of view, this contribution is very important, as it helps to prospect for hydrocarbon reservoirs, since δ¹¹B gives information on the organic matter maturation state of the oilfield.