Yedoma Permafrost Releases Organic Matter with Lesser Affinity for Cu2+ and Ni2+ as Compared to Peat from the Non-Permafrost Area: Risk of Rising Toxicity of Potentially Toxic Elements in the Arctic Ocean

Pollution of the Arctic Ocean by potentially toxic elements (PTEs) is a current environmental problem. Humic acids (HAs) play an important role in the regulation of PTE mobility in soil and water. The permafrost thaw releases ancient organic matter (OM) with a specific molecular composition into the Arctic watersheds. This could affect the mobility of PTEs in the region. In our study, we isolated HAs from two types of permafrost deposits: the Yedoma ice complex, which contains pristine buried OM, and the alas formed in the course of multiple thaw–refreezing cycles with the most altered OM. We also used peat from the non-permafrost region as the recent environmental endmember for the evolution of Arctic OM. The HAs were characterized using 13C NMR and elemental analysis. Adsorption experiments were conducted to assess the affinity of HAs for binding Cu2+ and Ni2+. It was found that Yedoma HAs were enriched with aliphatic and N-containing structures as compared to the much more aromatic and oxidized alas and peat HAs. The adsorption experiments have revealed that the peat and alas HAs have a higher affinity for binding both ions as compared to the Yedoma HAs. The obtained data suggest that a substantial release of the OM from the Yedoma deposits due to a rapid thaw of the permafrost might increase the mobility of PTEs and their toxicity in the Arctic Ocean because of much lesser “neutralization potential”.


Introduction
Environmental pollution by potentially toxic elements (PTEs) is nowadays an urgent global problem. PTEs are persistent environmental pollutants that do not degrade over time and are permanently present in the environment. Economic growth, together with transition to the green and renewable technologies and energy sources, has triggered a rapid increase in the consumption and production of copper, nickel, steel and other metals [1][2][3]. Emissions of copper and nickel smelter plants are significant sources of these elements to the environment [4][5][6]. This pollution is particularly critical in the sensitive Arctic ecosystems [7].
Although the sources of anthropogenic emissions are usually located outside the Arctic territories and their presence is caused by atmospheric transport from southern latitudes, significant industrial emissions from Cu and Ni smelters are observed in the Russian Arctic [8][9][10][11]. The levels of emission of inorganic elements into the atmosphere from the mining and processing plants of the Kola and Taimyr peninsulas of Russia were 350. 8 and 427 (Ni) and 798 and 474 (Cu) t/year, respectively, in 2001-2010 [12]. These amounts account for approximately 128% (nickel) and 196% (copper) of the annual emissions of these elements in Central Europe. Such high levels of PTE emissions have a negative impact on the flora and fauna of the region. While the highest concentrations of PTEs are found in the vicinity of the emitting sources and accumulate in bottom sediments, peatlands and permafrost soils, their mobilization poses a potential threat to the global environment [7,[13][14][15].
The mobility of metals in soil and water is regulated by hydrological parameters and soil regime (e.g., salinity, pH, temperature), redox conditions, chemical forms of the elements and organic matter (OM) [14,[16][17][18][19]. Among the many factors that affect PTE mobility, the complexation of metals with organic matter and specifically with humic substances (HSs) is one of the most crucial [16,20,21]. Humic substances are complex molecular systems, which are ubiquitous in the environment. HSs are defined based on their solubility characteristics. Humic acids (HAs) are not soluble at pH below 2, whereas fulvic acids (FAs) are soluble in the entire pH range [22]. PTEs form stable insoluble complexes with HAs by interacting with carboxylic, phenolic and heteroatom-containing functional groups of HAs [23][24][25][26]. The molecular composition of HAs varies significantly and depends on their origin. Many studies showed that molecular composition significantly affects the sorption capacity and stability of HA complexes with PTEs [27][28][29][30][31].
Peat and permafrost OM constitute a significant part of Arctic organic matter. Permafrost alone constitutes more than half of the Earth's soil organic carbon [32,33]. The molecular composition of peat unaffected by permafrost and permafrost OM can exhibit significant variations primarily due to the persistent low temperatures that restrict microbial transformations of permafrost OM [34][35][36]. However, recent climate warming triggered the mobilization of old organic carbon which was stored in permafrost for millennia. The rates of permafrost thawing are estimated on the order of 0.8-1.1 PgC per year [37]. This substantial release of organic carbon characterized by its unique molecular composition is now increasingly intruding into the Arctic watersheds, resulting in a shift towards older permafrost dominant OM [38,39]. This change in the molecular composition of OM could potentially lead to a change in PTE complexation and affect PTE presence and mobility in Arctic water and soil ecosystems. However, the knowledge regarding the complexation of permafrost humic acids with PTEs in the Russian Arctic, such as copper (Cu) and nickel (Ni), is very limited. Furthermore, the differences in the accumulation of PTEs compared to peat humic acids are not well understood.
Our study aims to fill the gap in the literature by investigating the impact of permafrost thawing on the mobility and complexation of potentially toxic elements with HAs in the Russian Arctic. We compare the complexation of Cu and Ni with humic acids derived from permafrost and peat soils. Our study focuses on a significant and novel factor: the shifting molecular composition of organic matter in Arctic watersheds resulting from permafrost thawing. This research aims to shed light on the potential impacts of permafrost thaw on the mobility and complexation of potentially toxic elements in Arctic ecosystems.

Sampling Site
The experiments were conducted on peat and permafrost HAs. A peat sample (PE) was collected from the ridge-lake complex of high moor Ilas sphagnum peat bog in the Arkhangelsk region of Russia (64.334170 • N, 40.609598 • E) from the 10-50 cm layer using a stainless-steel peat corer. Two permafrost samples (DY1 and DY2) were collected from the Duvanny Yar deposit near Cherskii settlement, the Republic of Sakha (Yakutia), Russia, using a UKB-12/25 rotary drill rig. The DY1 sample was collected from thermokarst depressions considered to be alas deposits characterized by multiple thaw-refreezing cycles (68. sample was collected and placed in a zip-lock plastic bag. The samples were frozen after collection and transported to the laboratory while still frozen. They were then stored at −25 • C in a freezer until lyophilization. In total, three samples, namely one peat sample (PE) and two permafrost samples (DY1 and DY2), were used in this study.
The characteristics of peat and permafrost samples such as pH, organic carbon (OC) and total nitrogen (TN) concentrations and initial moisture content (MC) are presented in Table 1.

Extraction of Humic Acids
The extraction of HAs was carried out from lyophilized peat and permafrost samples. Extraction of nonpolar organic components was conducted using a 1:1 v/v (hexane:ethanol) mixture in a Soxhlet extractor for 4 h. After that, the samples were dried in a vacuum oven at 40 • C for 12 h. Then, decalcification was carried out with 1 M HCl (1:10 solid:liquid solution phase weight ratio) until the pH value dropped down to 2. The acidified suspension was shaken for 6 h and left for 24 h [40]. The residual decalcified samples were neutralized by adding 1 M NaOH to pH 7, and then 0.1 M NaOH was added (1:10 sample:solution weight ratio) according to the International Humic Substances Society (IHSS) recommendations [41]. Then, NaCl was added until the final concentration of 0.3 M NaCl to the combined extract, which was centrifuged for 20 min at 10,000 rpm to coagulate colloidal particles. The supernatant was then separated from the precipitate and acidified to pH 1-2 with 1 M HCl. The samples were then left to precipitate for 24 h and centrifuged at 10,000 rpm. Purification of precipitated HA gels from silicates and ash components was carried out with a mixture of 0.1 M HCl and 0.3 M HF for 12 h. Then, the HA samples were purified by dialysis (3.5 kDa membrane) until a negative test for Cl − with AgNO 3 . The purified HA samples were freeze-dried, homogenized with mortar and pestle, and sieved through a sieve with d = 0.1 mm. All extraction and purification procedures, excluding dialysis, were performed under argon to prevent oxidation.

Characterization of the Initial Samples and Isolated Humic Acids
Quantitative 13 C solution state NMR spectra were acquired using an Avance-400 BioSpin spectrometer (Bruker, Bremen, Germany) operating at 400 MHz proton frequency. For the NMR study, 40 mg of HA was dissolved in 1 mL of NaOD in 40 wt. % in D 2 O, 99.5 atom % D (Sigma Aldrich, Germany), and 13 C NMR spectra were recorded according to the conditions described in [42,43]. To exclude the nuclear Overhauser effect, the INVGATE pulse technique was used. For complete carbon nuclei relaxation, the pulse delay was set to 8 s, and the number of scans averaged was 4500. All NMR spectra were acquired after centrifugation using a 5 mm broadband probe.
Moisture content was determined using an MX-50 Moisture Analyzer (AND, Tokyo, Japan). Soil pH was analyzed using the 9045D EPA method [44] using a SevenCompact pH/Ion meter and glass electrode InLab Expert Pro (Mettler Toledo, Columbus, OH, USA). OC and TN concentrations as well as mass concentrations of CHNS and O were analyzed using a 2400 Series II CHNS/O elemental analyzer (PerkinElmer, Waltham, MA, USA) according to [45].

Adsorption Experiments
The sorption kinetic experiments were carried out at 25 ± 0.5 • C, pH 5 ± 0.03 and 0.01 M KNO 3 . The weights of 15,000 ± 0.010 mg HA were measured using a BM-22 analytical balance (AND, Tokyo, Japan) and placed into 15 mL centrifuge tubes (Sarstedt, Nümbrecht, Germany). Aliquots of 10 mL of solution containing 20 mg/L Cu 2+ or 15 mg/L Ni 2+ were added into the tubes, which were placed onto an overhead shaker at 75 rpm (multi-rotator RS-60 (Biosan, Riga, Latvia)). The sorption was measured at different time intervals (5,10,15,25,30,45,60 and 90 min). Adsorption isotherm experiments were carried out at 25 ± 1 • C, pH 5 ± 0.03 and 0.01 M KNO 3 . Then, 15 mg of HA and 10 mL solution containing 1 -100 mg/L Cu 2+ or 0.5-100 mg/L Ni 2+ were shaken for 8 h at 25 ± 1 • C at 75 rpm. The effect of pH was assessed by adjusting the solution pH within the range of 2-6.5 using the same conditions as the isotherm experiments. pH adjustment was conducted using 0.1 M solutions of KOH and HNO 3 which were added dropwise to each solution.
Measurements were conducted using an axial ICP-AES 720-ES spectrometer (Agilent Technologies, Santa Clara, CA, USA) equipped with a low-flow axial quartz torch. The instrument was calibrated within the range of 0.025-10 mg/L. It was equipped with a 2.4 mm inner diameter injector tube (Glass Expansion, Melbourne, Australia), a double-pass glass cyclonic spray chamber (Agilent Technologies) and a OneNeb nebulizer (Agilent Technologies, Santa Clara, CA, USA). A Sc solution (20 ppm) prepared from 1000 mg/L Sc standard solution (High Purity Standards, North Charleston, SC, USA) was introduced online together with a sample solution using Trident Internal Standard Kit (Glass Expansion) and used as an internal standard. Emission line 361.383 nm was used for scandium as an internal standard correction. A 1 ppm QC standard solution was prepared from QCS-27-250 (High Purity Standards) and used to control accuracy and calibration stability after each batch of adsorption experiments. A 5% HNO 3 solution was used to control the blank signal to assess possible memory effects. The measurement results obtained on different lines were averaged with differences between them of less than 5% and were used for final calculations. All glass and polypropylene vessels were thoroughly rinsed using concentrated nitric acid and deionized water before adsorption and analysis procedures. Levels of all test elements in the blank samples were below their limits of detection (LODs).
Conditions of ICP-AES measurements are shown in Table 2. All lines were measured simultaneously (a MultiCal mode). Linear and rational (quadratic) functions were used for the calibration. The position and baseline corrections were performed for all the peaks of distinctive elements using the tools of the spectrometer software. Deionized water (18.2 MΩ × cm from a Milli-Q Academic system, Millipore, Molsheim, France) was used to prepare all the solutions and for washing.
All initial concentration solutions were treated the same way as samples to assess the possible sorption on vial walls and errors during initial dilution. The obtained concentration values were used as initial concentrations (C 0 ).
The pH of the solutions was analyzed using SevenCompact pH/Ion meter and glass electrode InLab Expert Pro (Mettler Toledo, Columbus, OH, USA) calibrated using 1.68, 4.01, 6.85 and 9.18 buffer calibration solutions (Hanna Instruments, Woonsocket, RI, USA).
The amount (mg) of Me 2+ ions adsorbed at equilibrium per gram of HA was calculated using the following equation: where q e is the quantity of Me 2+ ions adsorbed at equilibrium (mg/g), C 0 is the initial concentration of Me 2+ ions (mg/L), C e is the equilibrium concentration of Me 2+ ions (mg/L), V is the volume of the metal solution (L) and m is the weight of HA (g). The adsorption degree was calculated according to Equation (2): where η is the adsorption degree (%). The kinetics of Cu 2+ and Ni 2+ sorption was examined under batch conditions for a period of 90 min, and the acquired kinetic data were analyzed using kinetic models, namely the Lagergren pseudo-first-order (PFO) kinetic model [46], the Ho and McKay, or pseudosecond-order (PSO), kinetic model [47,48], and the Weber-Morris (W-M) intraparticle diffusion kinetic model [49], as shown in Equations (3)-(5), respectively.
where t is the time; q e and q t are the amounts of adsorbed metal ions per mass of sorbent (mg/g) at equilibrium and at time t, respectively; k 1 is the pseudo-first-order rate constant (min −1 ); k 2 is the pseudo-second-order rate constant (g·mg −1 ·min −1 ); k id is the intraparticle diffusion rate constant (g·mg −1 ·min −1/2 ); and C id is the parameter reflecting an effect of the boundary layer thickness (mg·g −1 ). The adsorption isotherms were constructed by fitting the experimental data to nonlinear Langmuir [50], Freundlich [51], Dubinin-Radushkevich (D-R) (equation adapted for adsorption from aqueous solution) [52] and Sips [53] adsorption models for both Cu 2+ and Ni 2+ . To fit these models and to calculate the parameters of the adsorption, the corresponding Equations (6)-(9) were used: where q e is the equilibrium and q m is the maximum adsorption capacity of metal ions (mg/g); C e is the equilibrium concentration of metal ions (mg/L); K L is the Langmuir constant (L/g); K F is the Freundlich constant (L/g); n is the indicator of the surface heterogeneity degree; β is activity coefficient which represents the adsorption energy (mol 2 ·kJ −2 ), from which mean free energy was calculated as E = 1/β 2 (kJ/mol); C s is the solubility (mg/L); R is the molar gas constant (J/(mol·K)); T is the temperature (K); K S is the Sips constant related to the energy of adsorption (L/g); and n is the Sips model exponent, which describes the heterogeneity of the adsorption sites on the sorbent surface. Adsorption kinetic curves and isotherms were plotted using the Python library matplotlib and SciPy. The same open-source software was used to determine sorption parame-ters [54]. To investigate the effect of different error functions on the isotherms measured, six error functions, namely r 2 , HYBRID, mean percent standard deviation (MPSD), average relative error (ARE), sum of squared errors (SSE) and average absolute error (EABS), were evaluated. The respective methods are described elsewhere [55]. The isotherm parameters were obtained by minimizing the error functions within the range of the studied concentrations. The analysis was conducted using the Python libraries and code mentioned above.
All experiments were conducted in triplicates, and the average values are reported.

Characterization of the Samples of Humic Acids Isolated in this Study
The results of structural group analysis ( 13 C NMR) and elemental composition of the three isolated HA samples are summarized in Table 3. Table 3. The carbon distribution among the structural fragments of HA obtained by 13 C NMR signal integration (% of total carbon); calculated ΣOx and C ar /C alk coefficients and elemental ratios obtained from the results of elemental analysis. According to the 13 C NMR data obtained, the PE and permafrost samples differed substantially in the content of oxygen-containing carboxylic, carbonylic and phenolic groups (we designated them as ΣOx) which are mostly responsible for metal ion complexation [56,57]. The HAs isolated from alas deposits (DY1) has the highest ΣOx value of 26.7% C, which is 1.5% higher than PE (25.2% C) and more than 5% higher as compared to DY2 (21.6% C).
However, the contents of carbohydrate fragments in the region of 55-112 ppm for both alas and Yedoma permafrost HAs DY1 and DY2 are 36% and 35%, respectively, which are higher than that for PE HA (33%). This difference in carbohydrate content can be attributed to the process of biodegradation, which is a primary pathway for carbohydrate decomposition. The peat sampling area is located in a warmer climate than permafrost. Higher temperatures, absence of permanently frozen layers, greater variety and abundance of microbial community lead to a higher rate of carbohydrate decomposition [58] and, as a result, a lower amount of carbohydrates in PE HAs.
The lowest overall distribution of C among oxygen-containing functional groups (Σ of 47-112 and 146-220 ppm region) was found for DY2 HA (63%), while PE HA (65%) and DY1 (68%) HAs had a relatively higher portion of C bonded to oxygen. This is supported by the elemental analysis data, specifically the O/C ratio which is higher for DY1 and PE (0.48 and 0.46) as compared to DY2 (0.41).
The aromaticity of HA C ar /C Alk was calculated as a ratio of spectral intensity of C atoms of aliphatic structures within a range from 0 to 112 ppm to the sum of aromatic structure C signals within a range from 112 to 166 ppm [59]. DY1 and PE were characterized by the higher values of C Ar /C Alk of 0.61 and 0.57, respectively, as compared to the sample of DY2 HA (0.50). This could be explained by the deeper oxidation and biodegradation of easily degradable aliphatic fragments in this permafrost sample [58,60,61]. The H/C ratio, which is also an indicator of aliphaticity [59], is approximately 10% higher for DY2 than for DY1 and PE HA (1.18, 1.09 and 1.07 respectively), which confirms their more aromatic structure.
The C/N ratio of 13.9 is significantly lower for DY2 HA, while PE and DY1 have C/N ratios of 17.1 and 17.2, respectively. The lower values of the C/N ratio and the higher content of nitrogen in the Yedoma permafrost HA sample may be connected with lower rates of microbial nitrogen consumption [62] or the changes in the vegetation during the formation of permafrost layers, as sphagnum spp. which are presently common for bogs have a higher C/N ratio than flora which was common in fens [63].
The concentration of another heteroatom in HA molecular composition, sulfur (S), differs insignificantly between all samples, with an S/C ratio of 0.011 ± 0.002.
We found that the alas and peat HA samples have the highest values of aromaticity and oxygen-containing functional groups. This indicates a higher proportion of strong metal-binding functional groups in the molecular composition of these two samples. The Yedoma sample is characterized by the more reduced character. This may be related to its significantly lower exposure to oxidation and biodegradation, as the Yedoma (DY2) permafrost deposit was below the active layer border and was not affected by thawing, oxidation and enhanced biodegradation in comparison with alas (DY1) and peat (PE) humic acids.

Effect of pH on Cu 2+ and Ni 2+ Adsorption on Humic Acids
The negative impact on the environment from Cu and Ni smelters not only is connected to the release of toxic levels of PTEs, but also brings about the high emissions of SO 2 gas, which may cause a significant decrease in the pH of the affected water and soil ecosystems [64][65][66][67]. A decrease in pH may affect the mobility of the elements in the environment due to its leaching from inorganic and organic complexes [68]. The typical pH values of the Arctic rivers and lakes vary in a wide range from 5 to 9 [69], while in Arctic regions affected by Ni and Cu smelter emissions, typical pH values drop down to 5-7 [7], and sometimes to 4 [70]. This is why we studied the effect of pH on Cu 2+ and Ni 2+ adsorption degree on isolated HAs in the pH range of 2 to 6.5. Figure 1 shows the dependence of adsorbed ions of elements at equilibrium on the pH of the solution.
for DY1 and PE HA (1.18, 1.09 and 1.07 respectively), which confirms their more aromatic structure.
The C/N ratio of 13.9 is significantly lower for DY2 HA, while PE and DY1 have C/N ratios of 17.1 and 17.2, respectively. The lower values of the C/N ratio and the higher content of nitrogen in the Yedoma permafrost HA sample may be connected with lower rates of microbial nitrogen consumption [62] or the changes in the vegetation during the formation of permafrost layers, as sphagnum spp. which are presently common for bogs have a higher C/N ratio than flora which was common in fens [63].
The concentration of another heteroatom in HA molecular composition, sulfur (S), differs insignificantly between all samples, with an S/C ratio of 0.011 ± 0.002.
We found that the alas and peat HA samples have the highest values of aromaticity and oxygen-containing functional groups. This indicates a higher proportion of strong metal-binding functional groups in the molecular composition of these two samples. The Yedoma sample is characterized by the more reduced character. This may be related to its significantly lower exposure to oxidation and biodegradation, as the Yedoma (DY2) permafrost deposit was below the active layer border and was not affected by thawing, oxidation and enhanced biodegradation in comparison with alas (DY1) and peat (PE) humic acids.

Effect of pH on Cu 2+ and Ni 2+ Adsorption on Humic Acids
The negative impact on the environment from Cu and Ni smelters not only is connected to the release of toxic levels of PTEs, but also brings about the high emissions of SO2 gas, which may cause a significant decrease in the pH of the affected water and soil ecosystems [64][65][66][67]. A decrease in pH may affect the mobility of the elements in the environment due to its leaching from inorganic and organic complexes [68]. The typical pH values of the Arctic rivers and lakes vary in a wide range from 5 to 9 [69], while in Arctic regions affected by Ni and Cu smelter emissions, typical pH values drop down to 5-7 [7], and sometimes to 4 [70]. This is why we studied the effect of pH on Cu 2+ and Ni 2+ adsorption degree on isolated HAs in the pH range of 2 to 6.5. Figure 1 shows the dependence of adsorbed ions of elements at equilibrium on the pH of the solution. From Figure 1, it can be observed that the equilibrium concentration of both elements decreased as the pH decreased for all three humic acid (HA) samples. The Cu 2+ adsorption degree decreased from its maximum value of (16 ± 1) mg/g at pH 6.5 down to 12-14 mg/g within a pH range from 3 to 6. An abrupt decrease was observed at pH 2: the corresponding qe accounted for 3.5-4.5 mg/g. On the contrary, Ni 2+ adsorption did not depend From Figure 1, it can be observed that the equilibrium concentration of both elements decreased as the pH decreased for all three humic acid (HA) samples. The Cu 2+ adsorption degree decreased from its maximum value of (16 ± 1) mg/g at pH 6.5 down to 12-14 mg/g within a pH range from 3 to 6. An abrupt decrease was observed at pH 2: the corresponding q e accounted for 3.5-4.5 mg/g. On the contrary, Ni 2+ adsorption did not depend significantly on pH in the range of 4 to 6.5 (the q e value varied from 8.1 up to 9.3 mg/g for all HA samples). Still, at pH < 4, there was an abrupt reduction observed down to 1.4-1.6 mg/g.
The increase in Cu adsorption degree at pH 6.5 (Figure 1a) can be attributed to the transformation of Cu species in solution at pH > 5 from Cu 2+ to CuOH + and Cu(OH) 2 and its precipitation at pH 6.5 where Cu(OH) 2 formation is somewhat more pronounced [71,72] . In contrast, Ni was presented as Ni 2+ in the whole studied pH range [73], and its adsorption was not significantly affected in the pH range from 4 to 6.5 (Figure 1b). The decrease in the adsorption degree of HAs at pH < 3 for both ions may be associated with the competition between the protonation of carboxylic groups of HAs which gradually dissociates at pH 2.5-7 and the adsorption of metal ions [31]. The results showed that the difference in the pH between 4 and 6 which is typical for acidified lakes and rivers near Cu and Ni smelters does not significantly affect the adsorption degree of these ions. Therefore, a pH 5 was selected for further adsorption studies.

Adsorption Kinetics and Isotherms
To assess the rate of Cu 2+ and Ni 2+ adsorption on humic acids, kinetic experiments were conducted. Figures S1 and S2 show the relationships between the amount of Cu 2+ and Ni 2+ adsorbed on HA and time. The experimental results were fitted using the PFO, PSO and W-M adsorption kinetic models. The calculated model parameters are outlined in Table 4. The adsorption process of Cu 2+ and Ni 2+ proceeds at relatively high rates. Data show a Cu-rich equilibrium for 30-45 min. The calculated q e was 14.4, 15.4 and 14.1 for PE, DY1 and DY2, respectively. The results for model fitting showed that the pseudo-secondorder model better describes the kinetics of Cu 2+ adsorption onto all three HA sorbents. This indicates that the limiting stage of Cu ion adsorption onto HAs is chemosorption. The adsorption rate of Ni 2+ was found to be even faster. The equilibrium concentration q e = 8.6 and 9.2 mg/g was achieved faster for PE and DY1 HAs, respectively, already after 10 min of adsorption reaction. For DY2, the equilibrium was achieved at 20 min with q e = 7.8 mg/g. These high rates of Ni adsorption negatively influenced the results of the kinetic modeling as there were limited points for curve fitting ( Figure S2). However, the fitting of the models showed that the PFO kinetic model describes Ni adsorption in each HA with r 2 0.89 for all HAs. The PFO model assumes that the initial adsorption rate is fast and that the rate of adsorption decreases as the amount of adsorbate on the HA surface approaches equilibrium. The rapid adsorption rates observed in our study are consistent with the findings of Lodygin et al., who reported that the adsorption equilibrium of Cd 2+ and Zn 2+ ions on permafrost humic acids is reached within 30 min [74]. On the basis of the obtained kinetic results, it was also concluded that 8 h is sufficient for the complete establishment of sorption equilibrium for both ions.
An adsorption isotherm study was conducted to determine the maximum adsorption capacity (q m ) of HSs for a particular metal ion and the affinity of HSs for that metal ion.
The calculated parameters for each model are presented in Table 5. It can be noted that for both ions and each HA sample, the determination coefficient r 2 was ≥0.90. This indicates that each model could be used to explain the adsorption behaviors. To properly evaluate the accuracy of a model in fitting experimental data, HYBRID, MPSD, ARE, SSE and EABS error functions together with r 2 were used. According to the obtained results, the Sips and Langmuir models demonstrated higher r 2 and lower error functions compared to the other models. The results of error function calculations are presented in Table S1.
Based on these findings, Sips and Langmuir models were selected as suitable models for describing the adsorption process of Cu 2+ on HAs. Figure 2 illustrates the adsorption isotherms for Cu 2+ and three HAs, confirming that both models provide a good fit to the experimental data. that for both ions and each HA sample, the determination coefficient r 2 was ≥0.90. This indicates that each model could be used to explain the adsorption behaviors. To properly evaluate the accuracy of a model in fitting experimental data, HYBRID, MPSD, ARE, SSE and EABS error functions together with r 2 were used. According to the obtained results, the Sips and Langmuir models demonstrated higher r 2 and lower error functions compared to the other models. The results of error function calculations are presented in Table  S1. Based on these findings, Sips and Langmuir models were selected as suitable models for describing the adsorption process of Cu 2+ on HAs. Figure 2 illustrates the adsorption isotherms for Cu 2+ and three HAs, confirming that both models provide a good fit to the experimental data.  In the case of the Langmuir adsorption model, it can be stated that the surface of HAs contains a finite number of identical sites available for adsorption, and the adsorption of Cu 2+ ions does not affect one another. According to the Langmuir model, Cu 2+ ions form a monolayer coverage on the surface of HAs without any interaction between the adsorbed In the case of the Langmuir adsorption model, it can be stated that the surface of HAs contains a finite number of identical sites available for adsorption, and the adsorption of Cu 2+ ions does not affect one another. According to the Langmuir model, Cu 2+ ions form a monolayer coverage on the surface of HAs without any interaction between the adsorbed ions [75,76]. At the same time, the Sips isotherm model, which combines the Langmuir and Freundlich models, is known to provide a better description of heterogeneous surfaces. Specifically, at low sorbate concentrations, it reduces to a Freundlich isotherm, while at high sorbate concentrations, it predicts a monolayer adsorption capacity consistent with the Langmuir isotherm [77,78].
The calculated q m parameters for Cu adsorption using both models show minimal differences, with Sips q m fluctuating around 28 ± 1 mg/g. However, the Langmuir model yields a higher q m value for DY2 (31 mg/g), whereas the q m values for PE and DY1 are approximately 27 mg/g. The K L and K S values for PE and DY1 HAs are also quite similar, with 0.56 and 0.55 for PE and 0.58 and 0.52 (L/g) for DY1, respectively. On the other hand, DY2 exhibits significantly lower K L and K S values, at 0.30 and 0.29 L/g, respectively. The higher q m and lower K L and K S values observed for DY2 suggest that while it has a slightly higher capacity for Cu 2+ adsorption, it also displays a lower affinity for Cu ions. This implies that DY2 may form less stable complexes with Cu 2+ compared to DY1 and PE.
It is worth noting that the adsorption degree of Cu ions by all humic acids (HAs) in the ratios of dissolved organic carbon (DOC) to Cu~100-1000:1 was observed to be 94 ± 1% for PE, 92 ± 2% for DY1 and 85 ± 5% for DY2. These ratios are specific to Arctic lakes and rivers that have been exposed to Cu pollution. In these environments, the concentration of DOC typically ranges from 1 to 10 mg/L or approximately 2 to 20 mg/L of HAs, while the concentration of Cu is around 1-20 µg/L [7,79,80].
The aforementioned observations suggest that DY2 exhibits a slightly lower adsorption capacity for Cu compared to PE and DY1. These findings align with the results of elemental and NMR analyses, which indicated higher concentrations of oxygen, higher O/C ratios, greater aromaticity and a higher proportion of C atoms bonded to oxygen for PE and DY1. In contrast, DY2 is more aliphatic and possesses a lower content of oxygen-containing groups that are primarily responsible for metal adsorption and chelation. The adsorption isotherms of Ni 2+ ions exhibit notable differences from the Cu 2+ curves (Figure 3). The error analysis results (Table S2) indicate that the Sips and D-R models provide better predictions for the adsorption of Ni compared to the Langmuir and Freundlich models. These two models demonstrate higher r 2 values and are characterized by lower error function values. A good fit of the D-R model for Ni ions suggests a complex mechanism of adsorption that involves both chemosorption and physical sorption in the pores of HAs. sistent with the Langmuir isotherm [77,78].
The calculated qm parameters for Cu adsorption using both models show minimal differences, with Sips qm fluctuating around 28 ± 1 mg/g. However, the Langmuir model yields a higher qm value for DY2 (31 mg/g), whereas the qm values for PE and DY1 are approximately 27 mg/g. The KL and KS values for PE and DY1 HAs are also quite similar, with 0.56 and 0.55 for PE and 0.58 and 0.52 (L/g) for DY1, respectively. On the other hand, DY2 exhibits significantly lower KL and KS values, at 0.30 and 0.29 L/g, respectively. The higher qm and lower KL and KS values observed for DY2 suggest that while it has a slightly higher capacity for Cu 2+ adsorption, it also displays a lower affinity for Cu ions. This implies that DY2 may form less stable complexes with Cu 2+ compared to DY1 and PE.
It is worth noting that the adsorption degree of Cu ions by all humic acids (HAs) in the ratios of dissolved organic carbon (DOC) to Cu ~100-1000:1 was observed to be 94 ± 1% for PE, 92 ± 2% for DY1 and 85 ± 5% for DY2. These ratios are specific to Arctic lakes and rivers that have been exposed to Cu pollution. In these environments, the concentration of DOC typically ranges from 1 to 10 mg/L or approximately 2 to 20 mg/L of HAs, while the concentration of Cu is around 1-20 µg/L [7,79,80].
The aforementioned observations suggest that DY2 exhibits a slightly lower adsorption capacity for Cu compared to PE and DY1. These findings align with the results of elemental and NMR analyses, which indicated higher concentrations of oxygen, higher O/C ratios, greater aromaticity and a higher proportion of C atoms bonded to oxygen for PE and DY1. In contrast, DY2 is more aliphatic and possesses a lower content of oxygencontaining groups that are primarily responsible for metal adsorption and chelation. The adsorption isotherms of Ni 2+ ions exhibit notable differences from the Cu 2+ curves ( Figure  3). The error analysis results (Table S2)  According to the results presented in Table 5, the maximum adsorption capacities obtained for Ni ions by D-R model follow the order PE (31) < DY1 (44) < DY2 (52) (mg/g). The calculated values of qe using Sips model are lower: 20, 32 and 27 mg/g for PE, DY1 According to the results presented in Table 5, the maximum adsorption capacities obtained for Ni ions by D-R model follow the order PE (31) < DY1 (44) < DY2 (52) (mg/g). The calculated values of q e using Sips model are lower: 20, 32 and 27 mg/g for PE, DY1 and DY2, respectively. The K S of Ni 2+ adsorption on PE HA is more than 2 times higher (0.41) than that for DY2 (0.19) and DY1 (0.17 L/g). This indicates that Ni has a higher affinity to PE than to the permafrost HAs. The mean free energy obtained by the calculation of D-R model parameters for both elements lies within 13-16 kJ/mol, which indicates that the adsorption of Ni and Cu on HAs follows a chemical ion-exchange mechanism [81].
The average values of Ni concentration in Arctic waterbodies affected by Ni pollution are generally lower compared to Cu in freshwater reservoirs of the Kola and Taymyr peninsulas; the average concentrations of Ni range from 0.9 to 4.3 µg/L, with levels reaching up to 147 µg/L near emission sources [7,80]. These concentrations correspond to ratios of approximately 100-1000:1 (OM:Ni). Under these conditions, the adsorption degree of HAs to Ni is 91 ± 5%, 89 ± 4% and 82 ± 3% for PE, DY1 and DY2, respectively, while at the ratio~10:1 at extremely low values of DOM concentration in water and/or high Ni pollution, the adsorption degree decreases to 21 ± 1% for all HA samples.
Overall, higher and more consistent q m values, as well as sorption constants, were observed for Cu compared to Ni when adsorbed onto the three HA samples using all four models employed in this study. Higher K S , K L and K F values were obtained for the PE HA sample, followed by DY1 and then by DY2. The latter sample has significantly lower metal binding affinity as compared to more oxidized OM.

Conclusions
We found that the HA samples isolated from the permafrost deposit had very different properties: the alas HA sample, which we defined as "the oxidized PM HA" (DY1), had properties very similar to the HA sample from peat, whereas Yedoma permafrost HA (DY2) was characterized by a much more reduced nature due to the higher contribution of the saturated aliphatic structures.
The adsorption experiments have revealed that the more oxidized and aromatic HA samples from the non-permafrost peat and from the alas permafrost had a higher affinity for binding both copper and nickel ions as compared to the well-preserved Yedoma permafrost HAs. The obtained data suggest that a substantial release of the organic matter from the Yedoma deposits due to a rapid thaw of the permafrost might bring about an increase in the mobility of PTEs and in their toxicity in the Arctic Ocean because of much lesser "neutralization potential". This factor should be taken into consideration for risk assessment scenarios while predicting the consequences of climate change on the pollution of the Arctic region.