domain: phys sde

text:
José Paulo Pinheiro

Elise Rotureau

Jerome F L Duval
email: jerome.duval@univ-lorraine.fr

Jérôme F L Duval

Addressing the electrostatic component of protons binding to aquatic nanoparticles beyond the Non-Ideal Competitive Adsorption (NICA)-Donnan level: Theory and application to analysis of proton titration data for humic matter

Keywords: Poisson-Boltzmann theory, NICA-Donnan model, Particle electrostatics, Proton titration curves, Nanoparticles, Humic matter

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Natural organic matter (NOM), the most important source of organic carbon on earth, is ubiquitous in soils and in aquatic ecosystems [1]. NOM refers to a broad spectrum of heterogeneous mixtures of macromolecules, nanoparticles and/or colloidal aggregates that differ in terms of physical and chemical properties including size, molecular weight, composition and charge [2]. The humic substances (HS), which represent the major NOM fraction, may be differentiated according to their solubility: the insoluble fraction, the humin, exists only in soils, whereas soluble fulvic and humic acids (hereafter denoted as FA and HA, respectively) are found in soils and natural waters [3]. Solubility of FA and HA further differs with respect to pH, with FA being fully soluble even under extreme acid conditions, and HA being soluble solely in weekly acidic and alkaline solutions [4]. Both FA and HA colloids are defined by nanometric dimensions with radius of ca. 1 nm and 2 to 10 nm [5], respectively. Accordingly, they are now generically termed as humic nanoparticles (HNPs in short) in literature.

The charge carried by HNPs strongly affects the mobility and the speciation of trace metal ions (M) in aquatic media. Indeed, it dramatically impacts on (i) the stability of HNPs-M complexes [6], [7]. (ii) the kinetics of the processes driving the formation and the dissociation of these complexes in solution [8], [9], and (iii) their lability, i.e. the extent to which complexes dissociate on the time scale of their diffusion toward M-consuming microorganisms [10], [11], [12]. The resulting connection between HNPs charge and metal bioavailability in aqueous environment motivates the need to address the electrostatic properties of HNPs at the proper quantitative level. The magnitude of the charge carried by HNPs inherently depends on solution pH and on background electrolyte concentration, both affecting the dissociation degree of the carboxylic and phenolic groups distributed throughout the HNPs volume [13]. In turn, the evaluation of HNPs charge properties from e.g. proton titration data analysis necessarily asks for an adequate differentiation between the chemical (intrinsic) and electrostatic components of the affinity of the protons for the HNPs carboxylic and phenolic sites [14]. Successful unravelling of these contributions is obviously tied to the robustness of the models in integrating the relevant electric double layer properties of HNPs, or, equivalently, the electrostatic potential distribution that is operational at the HNPs/solution interface.

So far, proton binding to HNPs has been mainly addressed from applications of the Windermere Humic Aqueous Model (WHAM) and of the Non-Ideal Competitive Adsorption (NICA) model developed in the 90's [15][16][17][18] . For the sake of conciseness, we refer the reader to Koopal et al. [START_REF] Koopal | Ion binding to natural organic matter: General considerations and the NICA-Donnan model[END_REF], [START_REF] Koopal | Equilibrium mono-and multicomponent adsorption models: From homogeneous ideal to heterogeneous non-ideal binding[END_REF] and Town et al. [7] for review and critical review of these modelling approaches. In their first attempt to integrate electrostatics within NICA formalism, de Wit et al. [START_REF] De Wit | Proton Binding to Humic Substances .1. Electrostatic Effects[END_REF] interpreted proton titration data measured on various humic substances at different solution ionic strengths by means of the Poisson-Boltzmann (PB) theory applied to hard particles, i.e. particles impermeable to ions from the background electrolyte.

Although such a PB-based modelling of humics electrostatic effects was reported to reproduce 'reasonably well' experimental titration data [START_REF] De Wit | Proton Binding to Humic Substances .1. Electrostatic Effects[END_REF], a simpler but empirical electrostatic representation was proposed where the potential is a priori considered constant within the whole particle body and in a surrounding solution volume adjusted to match proton titration data [13]. In its initial version, this 'Donnan-like' electrostatic model took the form of an equation involving 2 parameters describing empirically the dependence of a ˋDonnan volume´ on solution ionic strength [START_REF] Benedetti | Humic Substances Considered as a Heterogeneous Donnan Gel Phase[END_REF], replaced later by an equation involving a single adjustable parameter [15]. It is this last version of the ˋDonnan-like representation´ that is integrated in the NICA-Donnan (NICAD) model and, except for some differences detailed elsewhere [7], in the WHAM model.

Even though NICA-Donnan modelling is largely adopted in metal speciation and proton binding studies involving HNPs [16] or other types of nanoparticles [START_REF] Koopal | Ion binding to natural organic matter: General considerations and the NICA-Donnan model[END_REF], there are severe inconsistencies in its account of particle electrostatics, as recently detailed by e.g. Town et al. [7]. Briefly, the spatial distribution of the electrostatic potential at the HNPs/solution interface adopted in NICAD violates the fundamental Poisson-Boltzmann formalism, a funding element of the classical DLVO theory [START_REF] Lyklema | Fundamentals of Interface and Colloid Science Volume Iv Particulate Colloids Introduction to Colloid Science[END_REF]. Indeed, in the limit where particle radius well exceeds the electric Debye layer thickness in solution (a situation reached at sufficiently large electrolyte concentrations), the rigorous solution of the PB equation takes the form of a step-like potential profile where the potential is essentially constant inside the particle volume (this is the so-called Donnan potential) and zero outside: this is the true Donnan situation [START_REF] Duval | Electrokinetics of Diffuse Soft Interfaces. 1. Limit of Low Donnan Potentials[END_REF], [START_REF] Duval | Electrokinetics of Diffuse Soft Interfaces. 2. Analysis Based on the Nonlinear Poisson-Boltzmann Equation[END_REF], [START_REF] Overbeek | The Donnan Equilibrium[END_REF].

Accordingly, provided that medium salinity and particle size guarantee a correct application of the Donnan representation in line with the inequality κrp >>1 with rp the particle radius and 1/κ the Debye layer thickness, the Donnan volume cannot exceed the physical volume of the particles [7], [START_REF] Overbeek | The Donnan Equilibrium[END_REF]. This requirement is not respected by the NICAD model applied to HNPs [START_REF] Avena | Proton binding to humic acids: Electrostatic and intrinsic interactions[END_REF]. In addition, within the NICAD modelling framework, the Donnan electrostatic picture is adopted regardless of the size of the particles considered and of the background electrolyte concentration, i.e. irrespectively of the value taken by κrp.

This implies that the Donnan representation is basically considered valid over the whole spectrum of electric double layer regimes, from the thick to the thin double layer extremes corresponding to κrp <<1 and κrp >>1, respectively. Obviously, this assumption is incorrect, recalling that a Donnan potential profile is strictly operational under the only condition κrp >>1 [7]. Given a radius of say 5 nm for HA particles, it is estimated that κrp ∼ 1.7 to 5 in the range 10 mM to 100 mM electrolyte concentrations, respectively, this salinity range being commonly considered for HA proton titration experiments [START_REF] Milne | Generic NICA-Donnan model parameters for proton binding by humic substances[END_REF].

Therefore, under conditions where the electrostatic component is largest (i.e. at low κrp) or, equivalently, when screening of particle electrostatics is least significant, the condition underlying the applicability of a Donnan potential representation is not respected, which makes inappropriate the use of Donnan-based modelling framework for HNPs [7].

The above shortcomings of the Donnan electrostatic model applied to HNPs are not new. Already in 1999, Avena et al. [START_REF] Avena | Proton binding to humic acids: Electrostatic and intrinsic interactions[END_REF] explicitly recognized that "application of a simple Donnan model will only give results that are comparable to those obtained with the [PB] impermeable sphere model if the Donnan volumes used are larger than the specific volumes of the hydrated humics". In 2007, aware of the empirical nature of the adjustment of the 'particle Donnan volume' in NICAD, Companys et al. [START_REF] Companys | Electrostatic and specific binding to macromolecular ligands -A general analytical expression for the Donnan volume[END_REF] proposed a formulation of the Donnan volume on the basis of analytical solutions of Poisson-Boltzmann equation.

The approach, however, excludes salinity conditions of practical interest where Donnan representation for HNPs is not applicable and it does not integrate the whole PB-based spatial distribution of the electrostatic potential. In 2005, Duval et al. [START_REF] Duval | Humic substances are soft and permeable: Evidence from their electrophoretic mobilities[END_REF] reported an electrokinetic analysis of HNPs particles over a large range of pH and electrolyte concentration conditions on the basis of an electro-hydrodynamic theory developed for the generic case of soft particles (i.e. permeable to ions from background electrolyte) [START_REF] Duval | Electrophoresis of Diffuse Soft Particles[END_REF] whose electrostatic properties were tackled at the non-linear PB level. The authors thereby demonstrated the inapplicability of the hard particles-representation for HNPs that display, instead, the characteristic electrokinetic signature of soft particles whose intraparticulate potential approaches the Donnan potential at sufficiently large κrp. It is emphasized that the hard sphere-representation proposed earlier for HNPs [START_REF] De Wit | Proton Binding to Humic Substances .1. Electrostatic Effects[END_REF] is not compatible with the establishment of a Donnan potential in the particle body: this potential follows indeed from the compensation of the intraparticulate HNPs structural charges by electrolyte ions that have been transferred from bulk solution to particle volume by conduction-diffusion transport. Finally, Saito et al. [START_REF] Saito | Electrostatic potentials of humic acid: Fluorescence quenching measurements and comparison with model calculations[END_REF] revealed dramatic (sometimes exceeding 100 %) discrepancies between average HA particle potentials as determined by NICAD and that measured from fluorescence quenching. Instead, they evidenced better agreement with average potentials estimated on the basis of ion permeablesphere model and with adopting an electrostatic particle representation that integrates the extraparticulate electric double layer contribution.

In view of the above elements, a sound physicochemical modelling of HNPs protolytic properties where particle electrostatics is tackled at the proper non-linear Poisson-Boltzmann level for soft spheres is still lacking in literature. In this work, we report such a formalism and we elaborate a numerical procedure able to handle automatically, consistently and simultaneously, the recovery of proton titration curves measured on HNPs of given size at different solution ionic strengths. These results are systematically compared to those published using empirical NICAD model, thereby clearly identifying the cases where application of NICAD is erroneous or approximate. The numerical package designed in this work is available on demand.

MATERIALS

Polydisperse HNPs samples may form aggregates depending on pH, medium salinity and nature of the cations in solution [START_REF] Pinheiro | Dynamic properties of humic matter by dynamic light scattering and voltammetry[END_REF]. In turn, the evaluation of HNPs radius (rp) and molar mass (Mw) from proton titration data modelling remains uncertain as the reliability of the obtained results necessarily depends on the way particle electrostatics is pictured, recalling that rigorous PB-based potential distribution at HNPs/solution interface intrinsically depends on particle size. To provide a clear comparison -free of the above uncertainties-between HNPs-protons binding properties derived following our approach and those obtained from NICAD, the analysis below is based on protons titration data sets collected on HNPs for which particle radius is measured by independent techniques and for which NICAD-based fitting parameters are explicitly provided in literature. Using such particle size estimation, we evaluated Mw of the HNPs after exploiting Mw data reported for HNPs of composition similar than those considered in this work. As detailed below, this procedure constrains the modelling of protolytic titration data versus pH and electrolyte concentration and makes transparent the inherent electrostatic limits of NICAD model. In contrast, it advocates for the recourse to more rigorous soft PB theoretical framework. In line with the above elements, we selected five proton titration data sets published for HNPs with known particle radius, namely: the Ribeirão Preto Humic acid (FPHA) [START_REF] Botero | Isolation and purification treatments change the metal-binding properties of humic acids: effect of HF/HCl treatment[END_REF], Suwanee River Fulvic Acid (SRFA) [START_REF] Milne | Generic NICA-Donnan model parameters for proton binding by humic substances[END_REF], Laurentian Fulvic acid (LFA) [START_REF] Janot | PEST-ORCHESTRA, a tool for optimising advanced ion-binding model parameters: derivation of NICA-Donnan model parameters for humic substances reactivity[END_REF] and the fulvic (Hf:FA) and humic acid (Hf:HA) from Sierra de Buio, Spain [START_REF] Gondar | Effect of soil depth on acid properties of humic substances extracted from an ombrotrophic peat bog in northwest Spain[END_REF].

A detailed description of these HNPs, of their size measurement and of the estimation of their molar mass (Table 1) is provided in Supporting Material (SM, section A therein).

THEORY

The evaluation of HNP charge and protolytic properties is commonly carried out from the analysis of protons titration experiments performed in 1:1 electrolyte whose concentration is varied in the range 1 mM to 100 mM. The strategy thus consists in using the proton as a probe to explore the electrostaticdependent charging behaviour of HNPs. In turn, pending proper theoretical treatment, titration data provide the desired relevant electrostatic descriptors and charge of HNPs as a function of pH and electrolyte concentration.

Electrostatic and chemical descriptors of HNPs.

The total charge 𝑄𝑄 0 carried by a single HNP and accessible from titration experiments may be expressed as the integral over the intraparticulate volume of the density 𝜌𝜌 0 (𝑠𝑠) (in C m -3 ) of net charges stemming from the dissociation of the structural carboxylic and phenolic groups distributed throughout the HNPs volume:

𝑄𝑄 0 (pH, 𝑐𝑐 elec ∞ ) = 4𝜋𝜋 𝐹𝐹𝑉𝑉 𝑝𝑝 ∫ 𝜌𝜌 0 (𝑠𝑠) (𝑟𝑟, pH, 𝑐𝑐 elec ∞ )𝑟𝑟 2 d𝑟𝑟 𝑟𝑟 𝑝𝑝 0 (1)
where 𝑉𝑉 𝑝𝑝 = 4π𝑟𝑟 𝑝𝑝 3 /3 is the volume of an individual HNP (m 3 ), r is the radial coordinate with the particle centre positioned at r=0, 𝑟𝑟 𝑝𝑝 is the particle radius, F is the Faraday and 𝑄𝑄 0 is expressed as a concentration of monovalent charges per unit HNPs volume (mol m -3 ). Eq 1 makes explicit that 𝑄𝑄 0 and 𝜌𝜌 0 (𝑠𝑠) both depend on solution pH and on the background electrolyte concentration denoted as 𝑐𝑐 elec ∞ . The charge density 𝜌𝜌 0

(𝑠𝑠)
is determined by the local concentration of protons within the HNP body which itself depends on the potential profile 𝜓𝜓(𝑟𝑟). For the sake of convenience, we introduce the dimensionless local electrostatic potential defined by 𝑦𝑦(𝑟𝑟) = 𝑧𝑧𝐹𝐹𝜓𝜓(𝑟𝑟)/𝑅𝑅𝑅𝑅 where 𝑧𝑧 is the valence of the z:z background electrolyte (z=1 under commonly adopted titration experiment conditions). The radial dependence of 𝑦𝑦(𝑟𝑟) is governed by the Poisson-Boltzmann (PB) equation where contributions from both (immobile) structural charges carried by HNPs and (mobile) ions from the background electrolyte are accounted for. After straightforward arrangements, the dimensionless PB equation applicable to soft HNPs particles reads as [37] ‡ ‡ There is a typo error in eq 3 of ref. [START_REF] Rotureau | Relationship between swelling and the electrohydrodynamic properties of functionalized carboxymethyldextran macromolecules[END_REF] where a minus sign is missing in front of the charge density term.

𝑑𝑑 2 𝑦𝑦(𝑟𝑟) 𝑑𝑑𝑟𝑟̃2 + 2 𝑟𝑟̃𝑑𝑑 𝑦𝑦(𝑟𝑟) 𝑑𝑑𝑟𝑟̃- sinh[𝑦𝑦(𝑟𝑟)] = - 𝜌𝜌 0 (𝑠𝑠) �𝑟𝑟,pH,𝑐𝑐 elec ∞ � 2𝐹𝐹𝐹𝐹𝑐𝑐 elec ∞ ( 2
)
where 𝑟𝑟̃ is the dimensionless radial position defined by 𝑟𝑟̃= κ𝑟𝑟 with κ the reciprocal Debye length given by

𝜅𝜅 = �2𝐹𝐹 2 𝑧𝑧 2 𝑐𝑐 elec ∞ (𝑅𝑅𝑅𝑅𝜖𝜖 0 𝜖𝜖 𝑟𝑟 ) ⁄ (3)
with 𝜖𝜖 0 the dielectric permittivity of vacuum and 𝜖𝜖 𝑟𝑟 the relative dielectric permittivity of the medium. Eq 3 tacitly implies that the contribution of the protons and hydroxyl ions to the total solution ionic strength is not significant, which is legitimate under the pH and 1:1 electrolyte concentration conditions adopted in this work (and in many reports from literature [START_REF] Milne | Generic NICA-Donnan model parameters for proton binding by humic substances[END_REF]). Relaxing this condition is straightforward pending replacement (i) of eq 3 by

𝜅𝜅 = � ∑ 𝐹𝐹 2 𝑧𝑧 𝑗𝑗 2 𝑐𝑐 𝑗𝑗 ∞ (𝑅𝑅𝑅𝑅𝜖𝜖 0 𝜖𝜖 𝑟𝑟 ⁄ 𝑁𝑁 𝑗𝑗=1
) where the index j runs over all types of ions with bulk concentration 𝑐𝑐 𝑗𝑗 ∞ present in solution, including H + and OH -, and (ii) modification of the sinh term in eq 2 along the lines given elsewhere [START_REF] Zimmermann | Electrokinetics of a Poly(Nisopropylacrylamid-co-carboxyacrylamid) Soft Thin Film: Evidence of Diffuse Segment Distribution in the Swollen State[END_REF]. In addition, the current version of the program makes it possible the analysis of measurements performed in 2:1 electrolytes based on a proper transformation of eq 2 along the lines detailed in [6] (see eqs S1-S2 therein). Application of the formalism to other mixtures of electrolytes is further possible pending differentiation in eq 2 between the various Boltzmann terms associated with the respective contributions of the mobile ions to the overall charge of the extraand intra-particulate electric double layers. Last, the program allows, if required, the extraction of the Boltzmann factor for proton accumulation (or that of any ions present in the electrolyte) within the overall body of the particle. This is done by spatial integration of the relevant local Boltzmann factor, as detailed in [6]. Also, if required, the Boltzmann factor pertaining to ion accumulation in the extraparticulate double layer can be easily retrieved upon proper setting of the limits of the corresponding spatial integration.

The boundary conditions associated to eq 1 are provided by

y(𝑟𝑟̃→ ∞) = 0 ; 𝑑𝑑𝑦𝑦(𝑟𝑟) 𝑑𝑑𝑟𝑟̃� 𝑟𝑟=0 = 0 (4a,b)
where eq 4a reflects electroneutrality condition in bulk electrolyte solution and eq 4b results from symmetry of the potential profile at the particle centre. For the sake of demonstration, the developments below are illustrated for HNPs particle but the formalism remains valid for any particle type pending implementation of the relevant (Dirichlet-or Neumann-type of) electrostatic boundaries [START_REF] Duval | Dynamics of metal uptake by charged biointerphases: bioavailability and bulk depletion[END_REF] and consideration of the proper number of structural charges type (set to 2 for HNPs) [START_REF] Nederlof | Determination of Proton Affinity Distributions for Humic Substances[END_REF].

Under given pH and solution ionic strength conditions, the characteristic potential distribution 𝑦𝑦(𝑟𝑟) derived from eqs 1-2 takes the form of a bell shape with a maximum reached at the particle center, and it decreases gradually to zero value at few 1/κ distances far from the particle surface located at 𝑟𝑟 = 𝑟𝑟 𝑝𝑝 [START_REF] Rotureau | Relationship between swelling and the electrohydrodynamic properties of functionalized carboxymethyldextran macromolecules[END_REF].

In the limit κrp >>1, the electrostatic Donnan situation is achieved with 𝑦𝑦�𝑟𝑟 ≤ 𝑟𝑟 𝑝𝑝 � = 𝑦𝑦 D and 𝑦𝑦�𝑟𝑟 > 𝑟𝑟 𝑝𝑝 � = 0 where yD is the dimensionless Donnan potential [7]. It is emphasized that the local pH inside the HNPs volume depends on position r and thus differs from the solution pH to an extent that depends on the particle charge, or equivalently, on the magnitude of the intraparticulate potential.

The raw data 𝑄𝑄 obtained from proton titration experiments are classically expressed in equivalent concentration of monovalent charges per mass of titrated HNPs (i.e. Q is in mol kg -1 ). 𝑄𝑄 is then related to the volume charge density 𝑄𝑄 0 (mol m -3 ) of a single HNP particle (eq 1) according to 𝑄𝑄 = 𝑄𝑄 0 𝑉𝑉 𝑝𝑝 𝑁𝑁 𝐴𝐴 /𝑀𝑀 𝑊𝑊 = 𝑄𝑄 0 𝑣𝑣 𝑝𝑝 (5) where 𝑁𝑁 𝐴𝐴 is the Avogadro number, 𝑀𝑀 𝑊𝑊 (kg mol -1 ) is the molar mass of the considered HNP material, and 𝑣𝑣 𝑝𝑝 = 𝑉𝑉 𝑝𝑝 𝑁𝑁 𝐴𝐴 /𝑀𝑀 𝑊𝑊 is the specific volume of HNP (in m 3 kg -1 ). The conversion of 𝑄𝑄 0 to a charge expressed in standard experimental unit (mol kg -1 ) therefore requires knowledge of the particle radius rp and of the molecular weight 𝑀𝑀 𝑊𝑊 , or of the only specific volume 𝑣𝑣 𝑝𝑝 . This conversion is systematically required because PB-based electrostatics of particles does not depend on the total amount of charges carried by the particles but, instead, by their volume density of charges.

The structural charges carried by HNPs originate from the (de)protonation of carboxylic and phenolic groups. Due to the intrinsic chemical heterogeneity of HNP material, equilibrium protonation constants are classically modelled as continuous distributions [START_REF] De Wit | Proton Binding to Humic Substances .2. Chemical Heterogeneity and Adsorption Models[END_REF] involving mean protonation constants denoted as 𝐾𝐾 � 𝑎𝑎1 and 𝐾𝐾 � 𝑎𝑎2 for the carboxylic and phenolic groups, respectively, and the distribution width mH1 and mH2. The later reflects the chemical heterogeneity of the carboxylic and phenolic binding sites throughout HNPs due to differences in their molecular environments. The limits mH1=mH2=1 apply to ideal situations where HNP composition is homogeneous. The local structural charge density, 𝜌𝜌 0

(𝑠𝑠) (𝑟𝑟, pH, 𝑐𝑐 elec ∞ ) involved in eqs 1-2 is then provided by 𝜌𝜌 0 (𝑠𝑠) (𝑟𝑟, pH, 𝑐𝑐 elec ∞ ) = 𝜌𝜌 maxH1 (𝑠𝑠) 1+10 𝑚𝑚 H1 (p𝐾𝐾 � 𝑎𝑎1 -pH) exp[-𝑚𝑚 H1 𝑦𝑦(𝑟𝑟)] + 𝜌𝜌 maxH2 (𝑠𝑠) 1+10 𝑚𝑚 H2 (p𝐾𝐾 � 𝑎𝑎2 -pH) exp[-𝑚𝑚 H2 𝑦𝑦(𝑟𝑟)] (6)
where 𝜌𝜌 maxH1

(𝑠𝑠)
and 𝜌𝜌 maxH2

(𝑠𝑠)
are the maximum densities of structural charges (mol m -3 ) achieved under pH conditions where all carboxylic and phenolic sites are fully deprotonated, respectively. 𝜌𝜌 maxH1

(𝑠𝑠)
and 𝜌𝜌 maxH2

(𝑠𝑠) may be transformed in mol kg -1 , by multiplication by 𝑣𝑣 𝑝𝑝 , resulting in 𝑄𝑄 maxH1 and 𝑄𝑄 maxH2 , respectively. The set of eqs 1-6 fully determines the dependence of the searched charge 𝑄𝑄 on pH and ionic strength. A Fortran program (named SPBT for Soft Poisson Boltzmann-based Titration) was developed to compute eqs 1-6. In particular, the solution of eq 2 was derived according to the numerical strategy employed elsewhere [START_REF] Duval | Electrokinetics of Diffuse Soft Interfaces. IV. Analysis of Streaming Current Measurements at Thermoresponsive Thin Films[END_REF] using Fortran COLSYS subroutine [START_REF] Ascher | Collocation Software for Boundary-Value Odes[END_REF] which approximates the solution through spline-collocation at Gaussian nodes and selects mesh subdivision following an auto-adaptative strategy. This algorithm has already been proven extremely efficient for evaluation of e.g. complex soft particle electrokinetics [START_REF] Moussa | Remarkable Electrokinetic Features of Charge-Stratified Soft Nanoparticles: Mobility Reversal in Monovalent Aqueous Electrolyte[END_REF] and soft macrosurface electrohydrodynamics [START_REF] Duval | Electrokinetics of soft polymeric interphases with layered distribution of anionic and cationic charges[END_REF]. The one-component NICA equation commonly adopted for proton titration analysis is a Langmuir-Freundlich equation that involves, similarly to eq 6, the protolytic parameters 𝜌𝜌 maxH1

(𝑠𝑠)
, 𝜌𝜌 maxH2

Fitting proton titration data to theory: strategy.

The theoretical framework provided by eqs 1-6 involves 8 parameters that must be consistently adjusted to reproduce proton titration data collected at different solution ionic strengths, namely rp, Mw, 𝜌𝜌 maxH1

(𝑠𝑠)
, 𝜌𝜌 maxH2 (𝑠𝑠)

(or 𝑄𝑄 maxH1 and 𝑄𝑄 maxH2 ), 𝐾𝐾 � 𝑎𝑎1 , 𝐾𝐾 � 𝑎𝑎2 , mH1 and mH2, all of them being independent of the solution ionic strength provided that particles do not significantly swell with varying pH or 𝑐𝑐 elec ∞ , which is the case for e.g. FA and HA (see related discussion in the next section) [7], [START_REF] Hosse | Determination of electrophoretic mobilities and hydrodynamic radii of three humic substances as a function of pH and ionic strength[END_REF]. As the range of rp and Mw values is known (Table 1), this number basically reduces to 6, which is consistent with the conclusion by Lenoir and Manceau [START_REF] Lenoir | Number of Independent Parameters in the Potentiometric Titration of Humic Substances[END_REF] who evidenced, via

Principal Component Analysis (PCA) of proton titration data collected on a large series of fulvic and humic samples, that only 6 independent parameters could be adjusted to reproduce experiments with leastsquare minimization. The analysis of proton titration data using the current version of NICAD requires one additional parameter to be adjusted. This parameter, termed b, connects empirically the 'particle Donnan volume' VD (erroneously allowed to exceed Vp) where the potential is considered a priori constant, to the logarithm of the solution ionic strength, log I. Lenoir et al. [START_REF] Lenoir | Number of Independent Parameters in the Potentiometric Titration of Humic Substances[END_REF] argued that using NICAD, "the total covariance of the system cannot be reduced by measuring titration data at various ionic strengths because log I is correlated to log VD". This observation mirrors the failure of NICAD in establishing at the correct PB level the connection between potential profiles at the HNPs/solution interface and medium salinity. Janot et al. [START_REF] Janot | PEST-ORCHESTRA, a tool for optimising advanced ion-binding model parameters: derivation of NICA-Donnan model parameters for humic substances reactivity[END_REF] compared unconstrained and constrained proton procedures for titration data fitting with use of NICAD. They concluded that different sets of model parameters derived with different procedures could lead to similar quality of data fit, thereby questioning again the physical relevance of the obtained model parameters. Another important aspect is that the lack of obvious inflection points in the proton titration curves implies that the quality of the theoretical recovery of the titration data, regardless of the funding modelling approximations, may depend on the initial estimates of the adjusted parameters, which in turn conditions the performance of the algorithm to converge to a solution with minimisation of residuals [START_REF] Lenoir | Convergence-Optimized Procedure for Applying the NICA-Donnan Model to Potentiometric Titrations of Humic Substances[END_REF]. Accordingly, to obtain more reliable information on HNPs charge and electrostatic properties from confrontation between experiments and theory, the fitting procedure we adopt, detailed in the next section, considers simultaneously the reconstruction of the proton titration data (Q vs. pH) and that of its first derivative with respect to pH (dQ/dpH vs. pH), the so-called proton affinity spectrum.

Procedure for proton titration data reconstruction.

The parameters rp, Mw, 𝑄𝑄 maxH1 , 𝑄𝑄 maxH2 , 𝐾𝐾 � 𝑎𝑎1 , 𝐾𝐾 � 𝑎𝑎2 , mH1 and mH2 were optimized to reproduce the raw proton titration data and the associated affinity spectra by coupling our home-made Fortran program SPBT with PEST, an open source software for adjustment and estimation [49] of parameters by confrontation of theoretical models to experimental data.

The minimization procedure adopted in PEST for fitting theory to data is based on Gauss-Marquardt-Levenberg algorithm which minimizes the weighted sum Σ of squared differences between modelgenerated data (supplied here by SPBT) and corresponding experimental measurements. In the next section, we discuss the performance of SPBT-PEST in reproducing the Q and dQ/dpH vs.

pH data measured at different background electrolyte concentrations for the HNPs samples listed in Table 1. We further compare the corresponding sets of optimized parameters with those derived by NICADbased modelling and reported in literature for FPHA [START_REF] Botero | Isolation and purification treatments change the metal-binding properties of humic acids: effect of HF/HCl treatment[END_REF], SRFA [START_REF] Milne | Generic NICA-Donnan model parameters for proton binding by humic substances[END_REF], LFA [START_REF] Janot | PEST-ORCHESTRA, a tool for optimising advanced ion-binding model parameters: derivation of NICA-Donnan model parameters for humic substances reactivity[END_REF], Hf:FA [START_REF] Gondar | Effect of soil depth on acid properties of humic substances extracted from an ombrotrophic peat bog in northwest Spain[END_REF] and Hf:HA [START_REF] Gondar | Effect of soil depth on acid properties of humic substances extracted from an ombrotrophic peat bog in northwest Spain[END_REF] particles.

RESULTS AND DISCUSSION

Comparison between HNPs specific volume derived from SPBT-PEST and Donnan volume estimated from NICAD modelling

As detailed in the theoretical section, the particle radius rp and molecular weight 𝑀𝑀 𝑊𝑊 can be used to compute the particle specific volume 𝑣𝑣 𝑝𝑝 (m 3 kg -1 ). The Donnan volume VD significantly increases with decreasing solution ionic strength I. This dependence of VD on I, systematically observed with NICAD modelling regardless of sample types [START_REF] Koopal | Ion binding to natural organic matter: General considerations and the NICA-Donnan model[END_REF],

has been motivated by viscosity measurements and argued to stem from particle swelling [START_REF] Avena | Volume and structure of humic acids studied by viscometry: pH and electrolyte concentration effects[END_REF]. However, as discussed by Town et al., [7] this justification is questionable because it ignores the influence of primary and secondary electroviscous effects that originate from particle electric double layer interactions and interparticulate electrohydrodynamic interactions. These effects, in turn, can lead to substantial variations in the viscosity of particle dispersions with pH and ionic strength even under constant particle size conditions [START_REF] Zurita | The Primary Electroviscous Effect in Silica Suspensions -Ionic-Strength and Ph Effects[END_REF], [START_REF] Chan | Secondary Electroviscous Effect in a Charged Spherical Colloid[END_REF], [START_REF] Megias-Alguacil | The electroviscous effect in ethylcellulose latex suspensions. Effect of ionic strength and correlation between theory and experiments[END_REF]. Besides, may the particle size (and thus VD) increases significantly with decreasing I, one would anticipate a strong dependence of VD with pH, similarly to that observed for e.g.

responsive PNIPAM particles [START_REF] Martin | Structure of Multiresponsive Brush-Decorated Nanoparticles: A Combined Electrokinetic, DLS, and SANS Study[END_REF], essentially because the driving electrostatic force for particle swelling depends on particle charge and thus on both I and pH. In contradiction to this anticipation, VD is considered constant with pH in NICAD modelling. Last, for both HA and FA of soil and aquatic origins, particle size measurements performed with different techniques [7], [START_REF] Balnois | Atomic force microscopy of humic substances: Effects of pH and ionic strength[END_REF], [START_REF] Lead | Diffusion coefficients and polydispersities of the Suwannee River fulvic acid: Comparison of fluorescence correlation spectroscopy, pulsed-field gradient nuclear magnetic resonance, and flow field-flow fractionation[END_REF] show that there is no significant ionic strength-nor pH-dependent swelling/contraction of the particles in aqueous dispersions (the argument excludes extreme pH and I conditions where particles aggregation may take place). These elements highlight the empirical ad hoc fitting/calibration of VD for every salt concentration tested. In addition, for many of the situations considered in Table 1, VD exceeds 𝑣𝑣 𝑝𝑝 . This basically comes to position particle structural charges within a volume that exceeds the physical volume of the particle, which is inconsistent from an electrostatic point of view. Within the PB framework adopted in SPBT-PEST analysis of proton titration data and proton affinity spectrum, the HNPs specific volume 𝑣𝑣 𝑝𝑝 is independent of solution pH and ionic strength and, accordingly, it does not suffer from the inconsistencies of VD estimates listed above. Last, it is clear that for most of the salinity and particle size conditions examined, the κrp >>1 condition required for rigorous application of a Donnan representation of the electrostatic potential distribution at HNPs/solution interface is not satisfied (see Table 1), which further questions the physical meaning of the NICAD framework applied under such conditions. and the corresponding adopted model parameter values are listed in Table 2 together with their statistical dispersion evaluated by PEST. For the sake of conciseness, reconstruction of data pertaining to LFA, Hf:FA and Hf:HA are given in SM (section E, Figures S3-S5) and corresponding model parameters are provided in Table 2.

Performance of SPBT-PEST in recovering proton titration and affinity spectrum data.

Figure 1A shows a good agreement between titration experiments and SPBT-based theory for FPHA particles. In particular, the analysis reproduces with accuracy the electrostatic component of the proton binding to the particles, which is materialized by the correct positioning of the theoretical curves with varying background electrolyte concentration. This successful match between theory and experiments does not suffer from ad hoc adjustment of particle size and it intrinsically integrates the way intraparticulate and extraparticulate electrostatic potential distributions are modified with varying particle size and electrolyte concentration, features that are not shared by the empirical NICAD modelling framework. It is further emphasized that data treatment on the basis of SPB-PEST is significantly more constrained than that performed with using NICAD as the former, unlike the latter, considers simultaneous analysis of proton titration and affinity spectrum. Figure 1B shows some discrepancies between SPBT-PEST outcome and proton affinity spectrum data especially at low and high pH values. In the low pH range, suspensions of HNPs are indeed more prone to aggregation due to significant lowering of the charge carried by the particles [START_REF] Palmer | Dynamic light scattering measurements of particle size development in aqueous humic materials[END_REF], while pH measurements by glass membrane electrodes become unstable at sufficiently high pH due to unbalanced competition between protons and sodium ions from background electrolyte solution [START_REF] Covington | Procedures for Testing Ph Responsive Glass Electrodes at 25, 37, 65 and 85-Degrees-C and Determination of Alkaline Errors up to 1 Mol Dm-3 Na+, K+, Li+[END_REF]. In addition to these unavoidable experimental difficulties, the estimation of dQ/dpH may be subject to significant error as tiny deviations of the dependence of Q on pH compared to a smooth trajectory, may lead to significant changes in estimated dQ/dpH.

The titration data collected from literature for SRFA particles relate to pH lower than 7.5 (Figure 2A).

Accordingly, the chemical parameter values retrieved for the phenolic sites should be considered with caution as their dissociation is not significant in the pH range examined, which is supported by the absence of a second peak in the affinity spectrum (Figure 2B). The experimental titration data at 1 and 10 mM NaNO 3 are anomalously superimposed for this HNPs sample type, which may be explained by an inaccurate positioning of the experimental curves by Milne et al. [START_REF] Milne | Generic NICA-Donnan model parameters for proton binding by humic substances[END_REF]. The SRFA titration curve is marked by a steep increase of the particle charge with increasing pH from 3 to 6, similarly to that observed for a weak acid: this feature is consistent with the low log 𝐾𝐾 � 𝑎𝑎1 value of 3.2 derived for the carboxylic groups.

The lack of data at higher pH values is reflected by the large errors made in the estimation of the optimised parameters pertaining to the phenolic sites, particularly for 𝑄𝑄 maxH2 and mH2 (Table 2). The variations of the proton affinity spectrum at low pH are better defined for the sub-nanometric SRFA particles than those measured for the larger FPHA particles (Figure 1B), which supports that particles aggregationfacilitated with increasing particle size [START_REF] Palmer | Dynamic light scattering measurements of particle size development in aqueous humic materials[END_REF]-likely affects data quality in this pH range. The experimental results collected for LFA particles are of good quality up to pH 8 (SM, Figure S3) but some suspicious changes in the measurements can be observed at larger pH, especially at 30 mM ionic strength, which is particularly evidenced by inspection of the proton affinity spectra. The SPBT-PEST reconstruction of the LFA proton titration data is satisfactory for all electrolyte concentrations conditions examined, which contrasts with the theoretical quality of the corresponding proton affinity spectra. As mentioned above, this is so because the computation of the derivative dQ/dpH is very sensitive to changes of Q with pH, even very small. The results derived for Hf:FA (SM, Figure S4) and Hf:HA (SM, Figure S5) are quite similar, with the noticeable result that chemical heterogeneity viewed by the prism of mH1 is modest for both samples (Table 2). A difference is that carboxylic sites fraction is more significant in Hf:FA and associated log 𝐾𝐾 � 𝑎𝑎1 is significantly lower than that for Hf:HA. Inspection of Figures S4-S5

reveals that SPBT-based theory slightly overestimates experiments at low electrolyte concentrations (10 mM and 5 mM). Whereas it is difficult to provide firm explanations of such differences for Hf:FA and Hf:HA samples, possible origins of the mismatch between experiments and theory (mismatch that is more pronounced with NICAD-based modelling, see SM, section F) can be advanced: the effectiveness of the protocol for separating the sub-nanometric fulvic and humic fractions in these samples, an approximate positioning of the experimental titration data at 5 and 10 mM, or the occurrence of physico-chemical reactions affecting the structure of the particles at sufficiently low medium salinity.

Comparison between protolytic chemical parameters of HNPs retrieved by SPBT-

PEST and NICAD

As discussed above, the Donnan-like representation considered within NICAD modelling does not provide a consistent picture of the soft HNPs electrostatic properties. It is then of importance to address how the inconsistency of this unconstrained, empirical modelling approach is reflected in the NICADderived protolytic chemical parameters. For that purpose, under given medium salinity conditions. Surprisingly, NICAD manages to provide a reasonable estimates of proton titration data for FPHA and LFA at low ionic strengths where Donnan criterion is not satisfied (see Table 1): it does so, however, at the cost of unrealistic adjustment of VD already invoked by Town et al. [7] and Lenoir et al. [START_REF] Lenoir | Convergence-Optimized Procedure for Applying the NICA-Donnan Model to Potentiometric Titrations of Humic Substances[END_REF] Further comparison between 𝐾𝐾 � 𝑎𝑎1,2 derived by NICAD and SPBT-PEST (Table 2) shows important differences, especially for carboxylic groups for which this difference may amount to more than an order magnitude. This difference is most important for Hf:FA and Hf:HA samples for which NICAD modelling is inappropriate (SM, Figures S9-S10), and it remains marginal for SRFA and LFA despite of the corresponding poor NICAD-based titration data reconstruction (SM, Figure S7)

and/or the inapplicability of the Donnan criterion over the whole range of solution ionic strength considered (Table 1). For Hf:FA, the log 𝐾𝐾 � 𝑎𝑎1 (2.29±0.01) obtained by NICAD is further anomalously low for carboxylic acid and it compares to the first proton dissociation constant of a stronger acid like phosphoric acid (pKa=2.12) [START_REF] Sillen | Stability constants of metal-ion complexes, 2nd edition[END_REF]. As a matter of fact, this suspicious log 𝐾𝐾 � 𝑎𝑎1 value just mirrors the attempt of NICAD to balance its physically erroneous description of particle electrostatics via artificial adjustment of VD. The reader is referred to Table 1 and Table 2 to appreciate the over-or under-estimation of the quantities 𝑄𝑄 maxH1 , 𝑄𝑄 maxH2 , mH1 and mH2 by NICAD as compared to those derived by SPBT-PEST. It is worth mentioning that huge differences of these quantities (some approaching or exceeding 100%, see Unsurprisingly, adopting eq S3 provides better reconstruction of proton titration data than with using eq S4 simply because the degree of parameters adjustment offered to the user is larger with the former equation.

CONCLUSIONS

Since the 90's when the first versions of the Non-Ideal Competitive Adsorption-Donnan model (NICAD) were developed [15][16][17][18], the community of environmental physical chemists has largely employed the NICAD model to estimate the chemical charge carried by aquatic nanoparticles from proton titration data measured at different salt concentrations [15][16][17][18][START_REF] De Wit | Proton Binding to Humic Substances .1. Electrostatic Effects[END_REF][START_REF] Benedetti | Humic Substances Considered as a Heterogeneous Donnan Gel Phase[END_REF][START_REF] Milne | Generic NICA-Donnan model parameters for proton binding by humic substances[END_REF][START_REF] Botero | Isolation and purification treatments change the metal-binding properties of humic acids: effect of HF/HCl treatment[END_REF][START_REF] Janot | PEST-ORCHESTRA, a tool for optimising advanced ion-binding model parameters: derivation of NICA-Donnan model parameters for humic substances reactivity[END_REF][START_REF] Gondar | Effect of soil depth on acid properties of humic substances extracted from an ombrotrophic peat bog in northwest Spain[END_REF][START_REF] Nederlof | Determination of Proton Affinity Distributions for Humic Substances[END_REF][START_REF] Lenoir | Convergence-Optimized Procedure for Applying the NICA-Donnan Model to Potentiometric Titrations of Humic Substances[END_REF][START_REF] Tan | Proton binding to soil humic and fulvic acids: Experiments and NICA-Donnan modeling[END_REF][START_REF] Kouhail | Formation of mixed Eu(III)-CO3-fulvic acid complex: Spectroscopic evidence and NICA-Donnan modeling[END_REF][START_REF] Lodeiro | Acid-base properties of dissolved organic matter extracted from the marine environment[END_REF]. Accurate knowledge of the charge of these particles is fundamental as e.g. it largely contributes to define the fate, speciation and bioavailability of charged contaminants like metals in the aqueous compartments of ecosystems. In particular, there is a large amount of literature data reporting analysis by NICA-Donnan or related models, of the protolytic properties of humic matter nanoparticles (HNPs) commonly viewed as model surrogates for organic matter particles in aquatic environments or soils [18,[START_REF] Koopal | Ion binding to natural organic matter: General considerations and the NICA-Donnan model[END_REF][START_REF] De Wit | Proton Binding to Humic Substances .1. Electrostatic Effects[END_REF][START_REF] Benedetti | Humic Substances Considered as a Heterogeneous Donnan Gel Phase[END_REF][START_REF] Milne | Generic NICA-Donnan model parameters for proton binding by humic substances[END_REF][START_REF] Botero | Isolation and purification treatments change the metal-binding properties of humic acids: effect of HF/HCl treatment[END_REF][START_REF] Janot | PEST-ORCHESTRA, a tool for optimising advanced ion-binding model parameters: derivation of NICA-Donnan model parameters for humic substances reactivity[END_REF][START_REF] Gondar | Effect of soil depth on acid properties of humic substances extracted from an ombrotrophic peat bog in northwest Spain[END_REF][START_REF] Lenoir | Number of Independent Parameters in the Potentiometric Titration of Humic Substances[END_REF][START_REF] Lenoir | Convergence-Optimized Procedure for Applying the NICA-Donnan Model to Potentiometric Titrations of Humic Substances[END_REF][START_REF] Avena | Volume and structure of humic acids studied by viscometry: pH and electrolyte concentration effects[END_REF][START_REF] Tan | Proton binding to soil humic and fulvic acids: Experiments and NICA-Donnan modeling[END_REF][START_REF] Kouhail | Formation of mixed Eu(III)-CO3-fulvic acid complex: Spectroscopic evidence and NICA-Donnan modeling[END_REF][START_REF] Lodeiro | Acid-base properties of dissolved organic matter extracted from the marine environment[END_REF]. For the sake of illustration, ecotoxicologists often employ HNPs as model systems to analyse the speciation of hazardous metals in dispersions of charged colloids and the implications thereof in terms of identification of the toxic metal forms towards living organisms (e.g. plants, bacteria, or algae) [START_REF] Tsiridis | Interactive toxic effects of heavy metals and humic acids on Vibrio fischeri[END_REF].

Despite of its popularity, NICA-Donnan suffers from severe approximations on the way particle electrostatics is implemented, recalling that estimation of the charge carried by any colloidal particles necessarily requires a proper formulation of the electrostatic component of the binding of the chargedetermining ions (protons) to the ionogenic sites distributed at the particle surface and/or in the particle body volume. Namely, NICAD assumes the applicability of Donnan electrostatic representation for charged (nano)particles regardless of the ratio between particle size and thickness of the operational electric Debye layer thickness [7]. This approximation necessarily fails for HNPs with radius 1 nm to 10 nm and practical salt concentrations in the range 1 mM to 100 mM. A consequence is that NICAD requires, for proton titration data fitting, the adjustment of a 'particle Donnan volume' that (i) exceeds the physical size of the particles and (ii) systematically increases with decreasing salt concentration, even in the absence of significant particle swelling [7,[START_REF] Milne | Generic NICA-Donnan model parameters for proton binding by humic substances[END_REF][START_REF] Botero | Isolation and purification treatments change the metal-binding properties of humic acids: effect of HF/HCl treatment[END_REF][START_REF] Janot | PEST-ORCHESTRA, a tool for optimising advanced ion-binding model parameters: derivation of NICA-Donnan model parameters for humic substances reactivity[END_REF][START_REF] Gondar | Effect of soil depth on acid properties of humic substances extracted from an ombrotrophic peat bog in northwest Spain[END_REF]. Items (i) and (ii) are inconsistent from an electrostatic point of view and simply reflect the inapplicable electrostatic description given in NICAD for particles whose size does not significantly exceed the Debye layer thickness [7]. In line with these observations, force is to recognize that a sound alternative to NICAD modelling for NPs charge evaluation is urgently needed.

Accordingly, in this work we elaborate a new approach for the analysis of NPs proton titration data where Poisson-Boltzmann theory for soft particle electrostatics is consistently implemented. The particle charge is evaluated from proper spatial integration of the local particle charge density defined by an isotherm reflecting the chemical binding of charge determining protons to ionogenic sites located at the particle surface and within the particle body volume, and corrected by local electrostatics retrieved from numerical solution to the nonlinear Poisson-Boltzmann equation valid for weakly to highly charged particles. Unlike NICAD, this framework is used to generate proton titration curves with a consistent account of particle size-dependent potential profile over the whole spectrum of Debye layer extension,

i.e. from the thin to the thick double layer regime. A complete numerical package (that we make available for the community) combining a home-made FORTRAN program (named SPBT) for PB-based computations of particle electrostatics and charge, and a free module (called PEST) for optimized adjustment of the relevant chemical and electrostatic model parameters, is developed. The combined SPBT-PEST programs allow for (i) the automated and constrained reconstruction of measured proton titration curves and associated proton affinity spectra collected at different salt concentrations, and (ii) the evaluation of the relevant particle electrostatic and protolytic features.

The whole methodology is successfully tested for a range of HNPs particles whose titration data and size measurements are available in literature, and the performance of SPBT-PEST duo is further explicitly compared with that of NICAD. In detail, unlike NICAD [START_REF] Gondar | Effect of soil depth on acid properties of humic substances extracted from an ombrotrophic peat bog in northwest Spain[END_REF], SPBT-PEST successfully reproduces the dependence of measured particle charge on pH and salt concentration for all HNPs tested, including those for which Donnan representation is not applicable. For HNPs whose size legitimates the applicability of Donnan potential profile, NICAD succeeds in reproducing the experimental data albeit at the cost of a 'particle Donnan volume' that is adjusted to values that exceed the physical volume of the particles, which is physically impossible [7,[START_REF] Milne | Generic NICA-Donnan model parameters for proton binding by humic substances[END_REF][START_REF] Botero | Isolation and purification treatments change the metal-binding properties of humic acids: effect of HF/HCl treatment[END_REF][START_REF] Janot | PEST-ORCHESTRA, a tool for optimising advanced ion-binding model parameters: derivation of NICA-Donnan model parameters for humic substances reactivity[END_REF]. In addition, SPBT-PEST properly recovers proton titration data for HNPs with maintaining a constant particle radius for all salt concentrations tested, which contrasts with the empirical adjustment of the NICAD 'particle Donnan volume' found to systemically increase with decreasing salt concentration [START_REF] Milne | Generic NICA-Donnan model parameters for proton binding by humic substances[END_REF][START_REF] Botero | Isolation and purification treatments change the metal-binding properties of humic acids: effect of HF/HCl treatment[END_REF][START_REF] Janot | PEST-ORCHESTRA, a tool for optimising advanced ion-binding model parameters: derivation of NICA-Donnan model parameters for humic substances reactivity[END_REF][START_REF] Gondar | Effect of soil depth on acid properties of humic substances extracted from an ombrotrophic peat bog in northwest Spain[END_REF]. We find that the chemical parameters pertaining to proton binding and proton affinity distribution, as retrieved from data modelling by SPBT-PEST and NICAD, may largely differ, especially so for the smallest HNPs where Donnan representation is not appropriate.

This finding simply highlights the empirical electrostatic foundation of NICAD model and implications thereof in terms of particle charge estimation. NICAD may succeed to reproduce, quantitatively, experimental proton titration data even for particles whose size forbids a priori its application. It does so thanks to the adjustment of the 'particle Donnan volume', which artificially counterbalances the electrostatic deficiency of the NICAD model to achieve data fitting. As a support of this element, it is found that both the Donnan volume and the key chemical parameters involved in the model can vary by 100% depending on the reported form of the empirical NICAD equation that defines the dependence of 'particle Donnan volume' on medium salt concentration. In contrast, SPBT-PEST does not require the adjustment of a 'particle Donnan volume' for fitting proton titration curves, recalling that this volume, when relevant, necessarily identifies with the physical volume of the particle and the electrostatic potential therein (the Donnan potential) is independent of position. The Donnan electrostatic potential representation at the particle/solution interface is nothing else than the solution of the PB equation in the limit of thin double layers. As such, SPBT-PEST bypasses the necessity by NICAD to recourse to a questionable over-parameterisation of particle electrostatics to reproduce experimental data. Unlike NICAD [START_REF] Benedetti | Humic Substances Considered as a Heterogeneous Donnan Gel Phase[END_REF][START_REF] Milne | Generic NICA-Donnan model parameters for proton binding by humic substances[END_REF][START_REF] Duval | Electrokinetics of Diffuse Soft Interfaces. 2. Analysis Based on the Nonlinear Poisson-Boltzmann Equation[END_REF], the PB-based alternative we propose integrates consistently the intimate connection between particle double layer potential, particle size and background electrolyte concentration.

We believe that the here-highlighted shortcomings of NICAD in its phenomenological representation of particle electrostatics should be of concern for the community as NICAD, together with WHAM that suffers from similar electrostatic approximations, constitutes a widely employed basis for computing e.g.

metal-to-humic thermodynamic complexation or addressing metal bioavailability in relation to toxicity effects assessment. Unlike NICAD and WHAM, the PB theoretical framework developed in this study can be extended to account explicitly for the presence of multivalent ions in solution [6], an option that is particularly interesting for analysing particle charging behaviour in natural waters rich in divalent cations.

In addition, depending on the nature of particles considered, inhomogeneous particle structure composition and gradients in dielectric permittivity between particle body and extraparticulate medium can be easily implemented within the PB equation along the lines described elsewhere [START_REF] Duval | Dynamics of metal uptake by charged biointerphases: bioavailability and bulk depletion[END_REF], [START_REF] Duval | Impact of Electrostatics on the Chemodynamics of Highly Charged Metal-Polymer Nanoparticle Complexes[END_REF]. Such particle structure and electrostatic refinements may be valuable provided that independent measurements targeting these aspects (e.g. diffusion neutron scattering and dielectric spectroscopy) are available to avoid the adjustment of an ill-defined set of parameters when reconstructing particle proton titration curves. This research did not receive any specific grant from funding agencies in the public, commercial, or notfor-profit sectors.

Figure captions.

Tables.

Table 1: Comparison between HNPs specific volumes 𝑣𝑣 𝑝𝑝 derived in this work and model-generated Donnan volume VD computed from the parameters reported in literature for NICAD modelling. Eq S3 and eq S4 given in SM specify the way VD was estimated, i.e. either via the two-parameter (α, β) dependent or the one-parameter (b) dependent NICAD empirical equations that relates VD to the logarithm of ionic strength. Initial estimates (superscript i in the A. Detailed description of the HNPs considered in this work, of their size measurement and of the estimation of their molar mass.

A1. FPHA particle type. FPHA sample was extracted from a peat sample collected in the Mogi river region of Ribeirão Preto, São Paulo State, Brazil [1], following the procedure recommended by IHSS for extraction and purification of organic matter [2]. Particle diffusion coefficient was measured by Scanning stripping chronopotentiometry (SSCP, see details below). The quadruplicate measurements yielded values of 5.2, 4.1, 5.9 and 4.7×10 -11 m 2 s -1 . Using Stokes-Einstein equation to convert diffusion coefficient into particle radius rp, we obtained rp = 4.5 nm and a 95% confidence interval of 3.1 to 5.9 nm. The FPHA sample was extracted from a tropical peat and measured particle size shows that it is a relatively large HNP material. In the literature, the closest sample we found is PPHA, for which molar masses of 23 kg mol -1 has been measured by equilibrium UV scanning ultracentrifugation [3]. Perminova et al. [4] reported number (Mn) and weight (Mw) average molecular weights of 6.4 and 19.2 kg mol -1 , respectively, by size exclusion chromatography on 8 peat humic substances samples.

Determination of FPHA diffusion coefficient by SSCP. Diffusion coefficient was determined using Cd(II) as a probe according to the procedure detailed by Pinheiro et al. [5] The experiments were carried out using an Ecochemie µAutolab III and a PGStat 12 in conjunction with a Metrohm 663 VA stand (Metrohm, Switzerland). The setup was controlled by the GPES 4.9 software from EcoChemie, the Netherlands. A three electrodes configuration was used, with a Hg thin film plated onto a rotating glassy carbon disk (2 mm diameter, Metrohm) serving as a working electrode, a GC rod counter electrode, and an Ag/AgCl reference electrode from World Precision Instruments DRIREF-5 (electrolyte leakage < 8×10 -4 µL h -1 ). SSCP experiments were carried out in 20 mL of 10 mM NaNO3 solution with a Cd(II) concentration in the range 6 to 10×10 -7 mol l -1 and a FPHA concentration in the range 8 to 10×10 -5 kg l -1 .

The adopted experimental conditions were: deposition time (td) 45 s, oxidizing current (Is) 2×10 -6 A applied until the potential reached -0.300 V with respect to Ag/AgCl reference electrode, rotation speed 1000 rpm. All solutions were purged for 15 min at the beginning of every experiment and for 20 s after each measurement (assisted by mechanical stirring of the rotating electrode). Measurements were performed for a range of deposition potentials, from the foot to the plateau of the SSCP wave, i.e. from -0.850 to -0.600 V for Cd(II). m 2 s -1 , which correspond to particle radius in the range 0.85 to 1.05 nm. De Wit et al. [7] reported a value of 0.6 nm, placing this material in the low end of the 1-2 nm range commonly reported for fulvic acids.

Molar masses of 0.829 kg mol -1 and 0.545 kg mol -1 were reported for SRFA by Aiken et al. [8] and De

log 𝑉𝑉 𝐷𝐷 = α + β log 𝐼𝐼 (S3)
where α and β are parameters adjusted to match proton titration data and they depend on the nature of the HNP considered. The authors proposed later a 'simplification' of eq S3 with the following one-parameter dependent correlation between VD and log I [15] log

𝑉𝑉 𝐷𝐷 = 𝑏𝑏(1 -log 𝐼𝐼) -1 (S4)
Regarding the chemical component of the NICA equation, [16] it takes the form of a continuous equation that describes the competitive binding of cations to the HNP:

𝑄𝑄 𝑖𝑖,𝑇𝑇 = � 𝑛𝑛 𝑖𝑖 𝑄𝑄 maxH1 𝑛𝑛 𝐻𝐻1 � �𝐾𝐾 � 𝑖𝑖1 𝑎𝑎 𝑖𝑖,𝐷𝐷 � 𝑛𝑛 𝑖𝑖1 ∑ �𝐾𝐾 � 𝑖𝑖1 𝑎𝑎 𝑖𝑖,𝐷𝐷 � 𝑛𝑛 𝑖𝑖1 𝑖𝑖 �∑ �𝐾𝐾 � 𝑖𝑖1 𝑎𝑎 𝑖𝑖,𝐷𝐷 � 𝑛𝑛 𝑖𝑖1 𝑖𝑖 � 𝑝𝑝 1 1+�∑ �𝐾𝐾 � 𝑖𝑖1 𝑎𝑎 𝑖𝑖,𝐷𝐷 � 𝑛𝑛 𝑖𝑖1 𝑖𝑖 � 𝑝𝑝 1 + � 𝑛𝑛 𝑖𝑖 𝑄𝑄 maxH2 𝑛𝑛 𝐻𝐻2 � �𝐾𝐾 � 𝑖𝑖2 𝑎𝑎 𝑖𝑖,𝐷𝐷 � 𝑛𝑛 𝑖𝑖2 ∑ �𝐾𝐾 � 𝑖𝑖2 𝑎𝑎 𝑖𝑖,𝐷𝐷 � 𝑛𝑛 𝑖𝑖2 𝑖𝑖 �∑ �𝐾𝐾 � 𝑖𝑖2 𝑎𝑎 𝑖𝑖,𝐷𝐷 � 𝑛𝑛 𝑖𝑖2 𝑖𝑖 � 𝑝𝑝 2 1+�∑ �𝐾𝐾 � 𝑖𝑖2 𝑎𝑎 𝑖𝑖,𝐷𝐷 � 𝑛𝑛 𝑖𝑖1 𝑖𝑖 � 𝑝𝑝 2 (S5)
where Qi,T is the total amount of species i bound to HNPs (mol kg -1 ), 𝐾𝐾 � 𝑖𝑖1 and 𝐾𝐾 � 𝑖𝑖2 are the corresponding complexation parameters, ni1 and ni2 are the stoichiometric parameters relative to the binding of species i to the carboxylic (index 1) and phenolic (index 2) groups, respectively, 𝑄𝑄 maxH1 and 𝑄𝑄 maxH2 refer to the concentrations of carboxylic and phenolic groups, and p1 and p2 correspond to the intrinsic heterogeneity parameters of the carboxylic and phenolic binding sites. For the situation met in protolytic titrations experiments performed in inert NaNO3, it can be legitimately assumed that the monovalent cations do not compete with the hydrogen ion [17], and eq S5 then reduces to

𝑄𝑄 𝐻𝐻,𝑇𝑇 = 𝑄𝑄 maxH1 �𝐾𝐾 � 𝐻𝐻1 𝑎𝑎 𝐻𝐻,𝐷𝐷 � 𝑝𝑝 1 𝑛𝑛 𝐻𝐻1 1+�𝐾𝐾 � 𝐻𝐻1 𝑎𝑎 𝐻𝐻,𝐷𝐷 � 𝑝𝑝 1 𝑛𝑛 𝐻𝐻1 + 𝑄𝑄 maxH2 �𝐾𝐾 � 𝐻𝐻2 𝑎𝑎 𝐻𝐻,𝐷𝐷 � 𝑝𝑝 2 𝑛𝑛 𝐻𝐻2 1+�𝐾𝐾 � 𝐻𝐻2 𝑎𝑎 𝐻𝐻,𝐷𝐷 � 𝑝𝑝 2 𝑛𝑛 𝐻𝐻2 (S6)
Eq S6 used to perform the reported NICAD analysis of the titration data collected for the samples of interest in this work is identical to eq 6 with the equivalences mH1=nH1p1 and mH2=nH2p2.

C. Details of here-developed SPBT-PEST treatment of proton titration data and associated proton affinity spectra.

This section describes in detail the optimization of the model parameters rp, Mw, 𝜌𝜌 maxH1 In the first step of the data treatment, an input data file containing the measured HNPs charge Q (in mol kg -1 ) vs. pH at given solution ionic strength I and the associated proton affinity spectrum (dQ/dpH vs. pH) is created. This spectrum was computed using Matlab via cubic spline smoothing of the experimental Q vs. pH data and subsequent derivation of the splined data with respect to pH. The second step consists in establishing initial estimates for all parameters. As far as rp and Mw are concerned, their respective estimates are listed in Table 1 whereas those for 𝑄𝑄 maxH1 , 𝑄𝑄 maxH2 , mH1, mH2, 𝐾𝐾 � 𝑎𝑎1 and 𝐾𝐾 � 𝑎𝑎2 evaluated from the proton titration curve Q vs. pH and the corresponding affinity spectrum dQ/dpH vs. pH measured at the largest electrolyte concentration (so as to minimize particle electrostatics contribution) along the lines detailed in SM (section D, Figure S2).

The back and forth communication between SPBT and PEST to achieve theoretical reconstruction of titration data takes place in a third step. First (step 3a in Figure S1), the parameter template file instructs PEST to fit rp and Mw with the parameters 𝑄𝑄 maxH1 , 𝑄𝑄 maxH2 , mH1, mH2, 𝐾𝐾 � 𝑎𝑎1 and 𝐾𝐾 � 𝑎𝑎2 fixed to their respective initial estimates. PEST then creates the input parameter file and it calls for execution of SPBT that subsequently returns an output file containing the Q and dQ/dpH versus pH and I data generated for the set of parameters values given in the parameter file. PEST then reads this output file, computes the weighted sum Σ of squared differences between model-generated values and experimental data.

Subsequently, PEST uses a Gauss-Marquardt-Levenberg algorithm to generate an updated set of values for the adjustable parameters (those provided in the parameter template file), and it informs SPBT to reevaluate Q and dQ/dpH versus pH and I with this new set of parameters. PEST then reads the corresponding output file and compares the newly computed Σ with that obtained in the previous iteration.

The SPBT-PEST programs duo iterates this procedure until the relative parameter variation is smaller than a prescribed value after execution of a number of iterations specified by the user. It is stressed in PEST manual [18] that a relative parameters variation lower than 1% after 3 successive iterations is a good compromise between accuracy and running time, and these criteria were therefore adopted in our procedure. Once this parameters optimization criterion is satisfied, the updated values of rp and Mw are given by PEST in a result file that further contains the Q and dQ/dpH vs. pH and I theoretically generated by SPBT. In step 3b of Figure S1, SPBT-PEST duo resumes the full procedure detailed in step 3a albeit with now considering rp and Mw as non-adjustable parameters with values given by those obtained at the end of step 3a, and with adjusting 𝑄𝑄 maxH1 , 𝑄𝑄 maxH2 , mH1, mH2, 𝐾𝐾 � 𝑎𝑎1 and 𝐾𝐾 � 𝑎𝑎2 to achieve Q and dQ/dpH vs. pH and I data fitting. This is done via modification of the parameter template file that indicates to PEST which parameters it should now adjust and which ones should be maintained constant all along the iterated fitting procedure. Further details on the above procedure with step-by-step treated examples are provided in the form of a SPBT-PEST fitting manual that is available on request together with the SPBT executable file.

D. Evaluation of the initial estimates of the chemical parameters involved in eq 6.

An example of the procedure followed to evaluate initial estimates of the chemical parameters involved in eq 6 is provided in Figure S2 for FPHA for which proton titration data were published by Botero et al. [1] The raw data are presented in terms of charge equivalents per kg of humic matter versus pH. Titration data collected at large background electrolyte concentration (100 mM NaNO3 in Figure S2) are selected to minimize as much as possible the contribution of electrostatics to proton binding and thus to provide best initial estimates of the only chemical terms involved in eq 6. Initial estimates of 𝑄𝑄 maxH1 , 𝑄𝑄 maxH2 . The pH at which all carboxylic acids are titrated corresponds roughly to the onset of a plateau establishment for the particle charge Q with increasing pH. The clear location of this plateau is here necessarily impaired by the ongoing dissociation of the phenolic groups. For the lack of better, initial estimate of 𝑄𝑄 maxH1 is therefore read at the pH value marking a clear transition between dissociation features of carboxylic and phenolic HNP sites, i.e. at the pH value where the affinity spectrum displays a minimum with increasing pH (Figure S2B). Then, estimate for 𝑄𝑄 maxH2 can be obtained by simple difference between the maximum charge 𝑄𝑄 max,total and 𝑄𝑄 maxH1 defined in Figure S2.

Initial estimates of p𝐾𝐾 � 𝑎𝑎1 and p𝐾𝐾 � 𝑎𝑎2 .The two maxima in the proton affinity (Figure S2B) spectrum provides the desired initial estimates of p𝐾𝐾 � 𝑎𝑎1 and p𝐾𝐾 � 𝑎𝑎2 . Initial estimates of mH1 and mH2. The parameters mH1 and mH2 correspond roughly to the half-height widths of the peaks detected for the carboxylic and phenolic groups in the proton affinity spectrum. A difficulty for their evaluation is that experimental affinity spectrum is generally incomplete in terms of the covered pH range (especially at low and high pH values), which explains why initial guessed values for mH1 and mH2 are commonly set to 0.5. For situations where optimized mH1 and mH2 significantly differ from 0.5, we restarted the overall SPBT-PEST fitting procedure with choosing as initial guessed values of mH1 and mH2 those obtained from the previous optimization scheme performed with adopting 0.5 as initial estimates.

E. Performance of SPBT-PEST in recovering HNPs proton titration and affinity spectra data.

The model parameters adopted for SPBT-PEST analysis of the charge and electrostatic descriptors of LFA, Hf:FA and Hf:HA are collected in Tables 1-2 and we provide below in Figure S3 (LFA),

F. Comparison between SPBT-PEST and NICAD reconstructions of proton titration data.

For the sake of comparison, Figures S6-S10

  maxH1 and 𝑄𝑄 maxH2 ), 𝐾𝐾 � 𝑎𝑎1 , 𝐾𝐾 � 𝑎𝑎2 , mH1 and mH2 (see details of NICAD modelling in SM, section B). The electrostatic component of this NICA equation is tackled via a Donnan-like representation of the potential profile inside the particle body and in its aqueous surrounding, as commented in the Introduction section and further detailed in SM, section B.

  Figure S1 in SM (section C) summarizes the key steps of the proton titration data treatment by the coupled SPBT-PEST codes, which includes the estimations of sound initial guessed values for the chemical parameters 𝑄𝑄 maxH1 , 𝑄𝑄 maxH2 , mH1, mH2, 𝐾𝐾 � 𝑎𝑎1 and 𝐾𝐾 � 𝑎𝑎2 (SM, Figure S2, section D) as well as refinement of measured rp and Mw.Optimization of rp and Mw at one side and 𝑄𝑄 maxH1 , 𝑄𝑄 maxH2 , mH1, mH2, 𝐾𝐾 � 𝑎𝑎1 and 𝐾𝐾 � 𝑎𝑎2 at the other is performed according to a two-step process detailed in FigureS1. The reader is referred to SM (section C) for exhaustive details of the procedure driving the iterative dialog between PEST (in charge of the generating parameters and confronting experimental prediction to experimental data) with SPBT that computes (Q vs. pH and I) and (dQ/dpH vs. pH and I) data for given set of parameters.

Figures 1 and 2

 2 Figures 1 and 2 report the SPBT-PEST theoretical analyses of the proton titration and affinity spectrum data measured for FPHA and SRFA at different NaNO3 electrolyte concentrations, respectively,

Figure 1 :

 1 Figure 1: Proton titration curves for FPHA (A) and corresponding proton affinity spectra (B) collected at three NaNO3 electrolyte concentrations: 10 mM (∆), 30 mM (◊) and 100 mM (ο). Symbols: measurements. Dotted lines: SPBT-PEST modelling results. Model parameters are listed in Tables 1-2. The charge is expressed in moles of equivalent charges per kg of FPHA material.

Table 2 :

 2 SPBT-PEST optimized values of the chemical parameters involved in the Langmuir-Freundlich equation (eq 6) for the various HNPs types considered in this work. Indicated dispersions of the parameters correspond to the 95% confidence interval provided by PEST. The parameters are retrieved from SPBT-PEST and NICAD modelling of the proton titration curves measured at different salt reconstructions and reported in Figures S6-S10 for the five HNPs samples of interest. SPBT-PEST and NICAD theoretical reconstructions of these proton titration data versus pH and salt concentration are also provided in Figures S6-S10.

Figure 2 : 2 .

 22 Figure 2: Proton titration curves for SRFA (A) and corresponding proton affinity spectra (B) collected at three NaNO3 electrolyte concentrations (indicated). Symbols: measurements. Dotted lines: SPBT-PEST modelling results. Model parameters are given in Tables 1-2. The charge is expressed in moles of equivalent charges per kg of SRFA material.

  Supporting material. A. Detailed description of the HNPs considered in this work, of their size measurement and of the estimation of their molar mass. B. Details of the NICA-Donnan modelling framework. C. Details of here-developed SPBT-PEST treatment of proton titration data and associated proton affinity spectra. D. Evaluation of the initial estimates of the chemical parameters involved in eq 6. E. Performance of SPBT-PEST in recovering HNPs proton titration and affinity spectra data. F. Comparison between SPBT-PEST and NICAD reconstructions of proton titration data. Additional details on the procedure outlined in E with step-by-step treated examples are provided in the form of a SPBT-PEST user manual that is available on request together with the SPBT executable file.

Figure 1 :

 1 Figure 1: Proton titration curves for FPHA (A) and corresponding proton affinity spectra (B) collected at three NaNO3 electrolyte concentrations: 10 mM (∆), 30 mM (◊) and 100 mM (ο). Symbols: measurements. Dotted lines: SPBT-PEST modelling results. Adopted model parameters are listed in Tables 1-2. The charge is expressed in moles of equivalent charges per kg of FPHA material.

Figure 2 : 2 .

 22 Figure 2: Proton titration curves for SRFA (A) and corresponding proton affinity spectra (B) collected at three NaNO3 electrolyte concentrations (indicated). Symbols: measurements. Dotted lines: SPBT-PEST modelling results. Adopted model parameters are given in Tables 1-2. The charge is expressed in moles of equivalent charges per kg of SRFA material.

A2.

  SRFA particle type. The Suwanee river fulvic acid (SRFA) is one of the International Humic Substances Society (IHSS) reference materials and probably the most studied HNP sample type. In 2001, Hosse and Wilkinson [6] measured SRFA diffusion coefficient by Fluorescence correlation spectroscopy (FCS) under different pH and ionic strength conditions and reported values between 2.1 and 2.6 × 10 -10

  maxH1 , 𝑄𝑄 maxH2 ), 𝐾𝐾 � 𝑎𝑎1 , 𝐾𝐾 � 𝑎𝑎2 , mH1 and mH2 to reproduce the raw proton titration data and the associated affinity spectra by coupling SPBT with PEST. The optimization procedure is further summarized in the flowchart displayed in FigureS1and briefly discussed below.

Figure S1 :

 S1 Figure S1:Flowchart detailing the procedure adopted for the reconstruction of proton titration data and proton binding affinity spectra, and for the evaluation of the chemical and electrostatic properties of HNPs. See details in the text.

Figure S2 :

 S2 Figure S2: Illustration of the procedure followed for evaluation of the initial estimates the chemical parameters involved in eq 6 from inspection of (A) proton titration data and (B) corresponding proton affinity spectrum. The data (symbols) refer to FPHA humic acids in 100 mM NaNO3.

Figure S5 :

 S5 Figure S5: Proton titration curves for Hf:HA (A) and corresponding proton affinity spectra (B) collected at three NaNO3 electrolyte concentrations (indicated) Symbols: measurements. Dotted lines: SPBT-PEST modelling results. Model parameters are listed in Tables 1-2. The charge is expressed in moles of equivalent charges per kg of Hf:HA material.

  collect the experimental Q vs. pH and I data for the different HNPs samples of interest together with their theoretical fitting realized by NICAD and SPBT-PEST.

Figure S6 :

 S6 Figure S6: Proton titration curves for FPHA collected at three NaNO3 electrolyte concentrations (indicated). Symbols: measurements. Dotted lines: SPBT-PEST modelling. Full lines: NICAD modelling. Model parameters are listed in Tables 1-2. The charge is expressed in moles of equivalent charges per kg of FPHA material.

Figure S7 :

 S7 Figure S7: Proton titration curves for SRFA collected at three NaNO3 electrolyte concentrations (indicated) Symbols: measurements. Dotted lines: SPBT-PEST modelling. Full lines: NICAD modelling. Model parameters are listed in Tables 1-2. The charge is expressed in moles of equivalent charges per kg of SRFA material.

Figure S8 : 1 - 2 .

 S812 Figure S8: Proton titration curves for LFA collected at three NaNO3 electrolyte concentrations (indicated). Symbols: measurements. Dotted lines: SPBT-PEST modelling. Full lines: NICAD modelling. Model parameters are listed in Tables 1-2. The charge is expressed in moles of equivalent charges per kg of LFA material.

Figure S9 :

 S9 Figure S9: Proton titration curves for Hf:FA collected at three NaNO3 electrolyte concentrations (indicated). Symbols: measurements. Dotted lines: SPBT-PEST modelling. Full lines: NICAD modelling. Model parameters are listed in Tables 1-2. The charge is expressed in moles of equivalent charges per kg of Hf:FA material.

Figure S10 :

 S10 Figure S10: Proton titration curves for Hf:HA collected at three NaNO3 electrolyte concentrations (indicated). Symbols: measurements. Dotted lines: SPBT-PEST modelling. Full lines: NICAD modelling. Model parameters are listed in Tables 1-2. The charge is expressed in moles of equivalent charges per kg of Hf:HA material.

Table 1

 1

collects rp, 𝑀𝑀 𝑊𝑊 and 𝑣𝑣 𝑝𝑝 obtained from SPBT-PEST analysis of the proton titration data measured at different solution ionic strengths for the various HNPs of interest. For the sake of comparison, Table

1

further reports the NICAD-based Donnan Volume VD (also expressed in m 3 kg -1 ) estimated for the various HNPs types considered on the basis of NICAD model parameters published in literature.

Table 1 :

 1 Comparison between HNPs specific volumes 𝑣𝑣 Numbers between brackets correspond to κrp values, recalling that Donnan potential-representation is strictly applicable for κrp >>1.

	HNPs				Donnan fitting	Model-generated Donnan
	type	rp	Mw	𝑣𝑣 𝑝𝑝	parameters	volume VD (m 3 kg -1 )
		(×10 -9 m)	(kg mol -1 )	(m 3 kg -1 )	α,β for eq. S3	100 mM 10 mM	1 mM
					or b for eq. S4
	FPHA	3.1 to 5.9* ,i	19.2 to 23 i		α=0.60±0.32	0.0074	0.0138	0.0257
	[34]	5.0 ii	23 ii	0.0137	β=-0.27±0.05	(5.20)	(1.64)	(0.52)
	SRFA	0.6 to 1.05 i	0.5 to 1.5 i		b=0.87	0.0055	0.0407	0.3020
	[28]	0.6 ii	0.545 ii	0.0010		(0.60)	(0.20)	(0.06)

𝑝𝑝 derived in this work and model-generated Donnan volume VD computed from the parameters reported in literature for NICAD modelling. Eq S3 and eq S4 given in SM specify the way VD was estimated, i.e. either via the two-parameter (α, β) dependent or the one-parameter (b) dependent NICAD empirical equations that relates VD to the logarithm of ionic strength. Initial estimates (superscript i in the Table

)

and refined values (superscript ii) of particle size rp and molecular weight Mw by SPBT-PEST are further indicated.

*Determined in this work by Scanned Stripping Chronopotentiometry (SSCP), see details in SM (section A, therein).

Table 2

 2 collect 𝑄𝑄 maxH1 , 𝑄𝑄 maxH2 , mH1, mH2, 𝐾𝐾 � 𝑎𝑎1 and 𝐾𝐾 � 𝑎𝑎2 derived by NICAD and provided in the literature for the samples of interest, to be compared with the values obtained by SPBT-PEST. In addition, for the sake of completeness we report in Figures

S6-S10 (SM, section F) the NICAD-based fitting of the proton titration curves for all HNPs listed in

Table

1

together with those obtained by SPBT-PEST. Briefly, data reconstruction quality by NICAD for the largest FPHA and LFA HNPs considered in this work is quite comparable for all examined ionic strengths to that achieved by SPBT-PEST, whereas it is significantly worse for all other HNPs defined by lower rp. This conforms qualitatively with the fact that the NICAD assumption of a constant intraparticulate potential is best verified for large particles, more prone to satisfy the criterion κrp >>1

Table 2

 2 for LFA) simply reflects, again, the phenomenological nature of the NICAD fitting exercise that is based on an electrostatic representation unsupported by the particle size-dependent intra-and extra-particulate spatial distributions governed by PB formalism for soft particles. Even with application of NICAD model under conditions where its founding electrostatic approximation is not verified, proton titration data reconstruction can be more or less 'satisfactory': this underpins the empirical nature of NICAD-based fitting of such data.

for LFA data) are obtained by NICAD depending on the nature of the VD-log I relationship adopted (either one or two-parameters dependent relation, see eqs S3 and S4 in SM, respectively). This observation, further supported by the large differences in VD (see Table

1

  Table) and refined values (superscript ii) of particle size rp and molecular weight Mw by SPBT-PEST are further indicated. Numbers between brackets correspond to κrp values, recalling that Donnan potential-representation is strictly applicable for κrp >>1.

	HNPs				Donnan fitting	Model-generated Donnan
	type	rp	Mw	𝑣𝑣 𝑝𝑝	parameters	volume VD (m 3 kg -1 )
		(×10 -9 m)	(kg mol -1 )	(m 3 kg -1 )	α,β for eq. S3	100 mM 10 mM	1 mM
					or b for eq. S4
	FPHA	3.1 to 5.9* ,i	19.2 to 23 i		α=0.60±0.32	0.0074	0.0138	0.0257
	[34]	5.0 ii	23 ii	0.0137	β=-0.27±0.05	(5.20)	(1.64)	(0.52)
	SRFA	0.6 to 1.05 i	0.5 to 1.5 i		b=0.87	0.0055	0.0407	0.3020
	[28]	0.6 ii	0.545 ii	0.0010		(0.60)	(0.20)	(0.06)
					b=0.35 ± 0.02	0.0004	0.0009	0.0019
	LFA	1.9 to 3.6 i	5.5 to 11 i			(3.74)	(1.18)	(0.37)
	[35]	3.6 ii	10 ii	0.0118	α=0.93±0.20	0.0158	0.0295	0.0550
					β=-0.27±0.05	(3.74)	(1.18)	(0.37)
	Hf:FA	0.75 to 1.0 i	0.5 to 1.5 i		b=0.43	0.0007	0.0019	0.0052
	[36]	1.0 ii	0.9 ii	0.0028		(1.04)	(0.33)	(0.10)
	Hf:HA	0.75 to 1.0 i	1.0 to 3.0 i		b=0.34	0.0005	0.0010	0.0023
	[36]	1.0 ii	1.8 ii	0.0014		(1.04)	(0.33)	(0.10)

*Determined in this work by Scanned Stripping Chronopotentiometry (SSCP), see details in SM (section A, therein).

Table 2 :

 2 SPBT-PEST optimized values of the chemical parameters involved in the Langmuir-Freundlich equation (eq 6) for the various HNPs types considered in this work. Indicated dispersions of the parameters correspond to the 95% confidence interval provided by PEST. The parameters are retrieved from SPBT-PEST and NICAD modelling of the proton titration curves measured at different salt reconstructions and reported in Figures S6-S10 for the five HNPs samples of interest. SPBT-PEST and NICAD theoretical reconstructions of these proton titration data versus pH and salt concentration are also provided in Figures

S6-S10.

	Hf:FA SPBT-PEST	4.35 ± 0.09	3.81 ± 0.01 0.91 ± 0.02 1.75 ± 0.06	7.1 ± 0.3	0.33 ± 0.01
		NICAD/eq.S4	5.0±0.01	2.29±0.01	0.44	0.94 ±0.03	6.85±0.02	0.59
		[36]
	Hf:HA SPBT-PEST	3.00 ± 0.05	4.40 ± 0.02 0.83 ± 0.03	1.3 ± 0.2	7.5 ± 0.2	0.49 ± 0.05
		NICAD/eq.S4	3.44± 0.02	3.35 ± 0.01	0.46	0.64 ± 0.04 6.68 ± 0.03	0.88
		[36]
	NA (not available)*: NICAD fit was carried out only for site 1 type (i.e. carboxylic) due to the experimental data
	covering a reduced pH range (see Fig. 2)
	HNPs	Model	𝑄𝑄 maxH1	log 𝐾𝐾 � 𝑎𝑎1	mH1	𝑄𝑄 maxH2	log 𝐾𝐾 � 𝑎𝑎2	mH2
	type		(mol kg -1 )			(mol kg -1 )
	FPHA SPBT-PEST	4.08 ± 0.02	4.26 ± 0.02 0.46 ± 0.01 1.39 ± 0.07 8.57 ± 0.09 0.59 ± 0.03
		NICAD/eq.S3	3.18 ± 0.92	3.65 ± 0.23	0.66	3.02 ± 2.4	8.07 ± 0.77	0.29
		[34]
	SRFA SPBT-PEST	5.47 ± 0.06	3.20 ± 0.02 0.52 ± 0.01	1.0 ± 0.4	7.8 ± 0.2	0.5 ± 0.2
		NICAD/eq.S4	5.55	3.01	0.42	NA*	NA*	NA*
		[28]
	LFA	SPBT-PEST	3.79 ± 0.12	3.94 ± 0.03 0.52 ± 0.02 4.94 ± 0.10	8.5 ± 0.2	0.27 ± 0.01
		NICAD/eq.S4	2.8 ± 1.2	2.5 ± 0.2	0.62 ± 0.17	7.4 ± 2.5	7.4 ± 0.7	0.18 ± 0.06
		[35]
		NICAD/eq.S3	5.1 ± 0.0	4.1 ± 0.1	0.41 ± 0.01	2.8 ± 0.0	8.6 ± 0.1	0.46 ± 0.02
		[35]

SUPPORTING MATERIAL

Addressing the electrostatic component of protons binding to aquatic nanoparticles beyond the Non-Ideal Competitive Adsorption (NICA)-Donnan level: theory and application to analysis of proton titration data for humic matter.

José Paulo Pinheiro, 1 Elise Rotureau, 1 Jérôme F. L. Duval 1, * 1 Université de Lorraine, CNRS, Laboratoire Interdisciplinaire des Environnements Continentaux (LIEC), UMR 7360, Vandoeuvre-lès-Nancy, F-54000, France.

* Corresponding author: jerome.duval@univ-lorraine.fr Tel: 00 33 3 72 74 47 20 This supporting material contains 12 pages, 6 equations, 18 references and 10 figures. It is organized as follows:

A. Detailed description of the HNPs considered in this work, of their size measurement and of the estimation of their molar mass.

B.

Details of the NICA-Donnan modelling framework.

C. Details of here-developed SPBT-PEST treatment of proton titration data and associated proton affinity spectra.

D.

Evaluation of the initial estimates of the chemical parameters involved in eq 6.

E.

Performance of SPBT-PEST in recovering HNPs proton titration and affinity spectra data.

F. Comparison between SPBT-PEST and NICAD reconstructions of proton titration data.

Wit et al. [7], respectively, while Her et al. [9] reported value of 1.116 kg mol -1 obtained by HPSEC/DOC and 1.385 kg mol -1 by HPSEC/UVA measurements.

A3. LFA particle type.

LFA is a soil fulvic acid, extracted and purified from a podzol obtained in the Laurentian forest preserve of Laval University, Quebec, Canada [10]. The diffusion coefficients measured by voltammetry for the LFA range from 6×10 -11 to 1.2×10 -10 m 2 s -1 at 10 mM ionic strength, corresponding to particle radius from 1.9 to 3.6 nm [11]. The LFA was extracted from a forest soil and its size is larger than that of SRFA. Shinozuka et al. [12] measured Mn of 9.7 and 11 kg mol -1 for Inogashira and Elliot soil fulvic acids by MALDI-TOF MS, while Perminova et al. [4] reported Mn of 5.5 kg mol -1 and Mw of 10.8 kg mol -1 averaged over 9 different samples.

A4. Hf:FA and Hf:HA particle types. Hf:FA and Hf:HA are the fulvic and humic acids extracted from the upper horizon (0-60 cm) of an ombrotrophic peat bog in Sierra del Buio, Northwest Spain. Pinheiro et al. [5] reported similar diffusion coefficients derived by FCS and SSCP, for both particle types with values in the range 2.3-2.8×10 -10 m 2 s -1 at 10 and 100 mM ionic strength, which corresponds to rp between 0.75 and 1.0 nm. This small particle size may be due to the fact that this material was extracted from the soil upper horizon where humification process is still recent. According to the authors [13], the collected sample contained reddish-brown fibric materials, while in the deeper horizons plant detritus was fully decomposed [13]. Hf:FA and Hf:HA particle radius indicates that they are a peculiar sample. Such a small radius is not compatible with large molar masses, thus we defined a low limit of 0.5 kg mol -1 for the FA and 1.0 kg mol -1 for the HA and fixed the higher limit at 1.5 and 3 kg mol -1 for the FA and HA samples, respectively in line with the Mw reported for small fulvic and humic NPs [7].

B. Details of the NICA-Donnan modelling framework.

NICA-Donnan model assumes a priori that the potential within the whole particle volume is constant (and takes the Donnan value ψD), and it further considers the accumulation of the only counterions within this volume and it neglects that of the coions that is however operational at sufficiently low particle charge and/or large electrolyte concentrations. In this model, cations accumulate inside the volume of (negatively charged) HNPs particles according to the Boltzmann accumulation factor

where ai and zi are the bulk activity and valence of cation i, respectively, F is the Faraday constant, R the gas constant and T is the temperature. Within the Donnan representation, the structural charge Q (in mol/kg of HNP material) located in the HNPs is neutralized by the cations, and the Donnan volume VD (which identifies with the specific hydrodynamic volume of the particle in m 3 /kg of HNP material) is defined by the electroneutrality equation:

Benedetti et al. [14] proposed an empirical relationship between VD and the solution ionic strength (I) for HNPs: Model parameters are listed in Tables 1-2. The charge is expressed in moles of equivalent charges per kg of Hf:FA material.