Conductometric evidence for accurate dissociation behaviour of surface groups on nanorods

Abstract

The surface charge of nanoparticles is usually quantified under the assumption that strong electrolytes fully dissociate in aqueous systems. Here, we demonstrate the effect of counterion condensation on charge through meticulous corrections to the individual dissociation constants. By considering the electrostatic potential energy of charged moieties on the particle surface and by applying Poisson–Boltzmann statistics to assess the electrostatic double layer surrounding the particles, we were able to simulate experimental potentiometric and conductometric titration curves for both solutions and colloidal systems with decent correlation. As such, we can predict that due to counterion condensation, our chosen rod-like particles (phosphorylated cellulose nanocrystals) never carry more than 0.5 mmol g⁻¹ charge, amounting to a mere fraction of the theoretical maximum of 2 mmol g⁻¹. In addition, we gained insight into the role of carbonate impurities in the system, which reliably indicate the first equivalence point of the highly charged nanorods by causing the brief rise in conductivity before a plateau is reached. We therefore discourage the interpretation of conductometry curves for nanoparticles by traditional linearization. Instead, we propose a more reliable method for nanoparticle charge determination from the conductometric data based on calculating the fraction of counterions entering the bulk solution instead of merely observing conductivity or pH.

---

Introduction

Nanoparticle charge in aqueous systems determines its colloidal stability, its capacity for ion exchange, its behavior in an electric field, and many of its self-assembling properties.¹ Generally, the surface charge is generated from adsorbed or covalently bound electrolyte moieties, which are prone to dissociation in aqueous systems.^2,3 If these moieties are confined in close proximity to each other on a particle surface, however (compare Fig. 1b as opposed to Fig. 1a), the electrostatic potential of the surface charges can prevent further dissociation.^4–7 In this case, their overlapping electrostatic potentials jointly impede the further abstraction of counterions and thereby shift the dissociation equilibrium of the chargeable moieties.^5,8 A qualitative comparison of the free energies of the dissociation reactions for ions and surface groups is shown in Fig. 1.

Fig. 1 While dissolved ions dissociate fully (a), surface charges do not (b), as their proximity gives rise to enhanced electrostatic interactions, turning the exergonic dissociation of ions (c) into a more endergonic process (d).

For polymeric polyelectrolytes, which can be envisioned as charged lines, this phenomenon of suppressed dissociation has been mathematically described by Manning and termed counterion condensation.⁹ Independent of the actual dissociating species, Manning highlighted the distance between two distinct charges along the linear polyelectrolyte as crucial, with sufficient overlap of the individual potentials giving rise to counterion condensation.^9–11

This general concept applies to all dissociation equilibria contributing to the net surface charge, whether the charges are lined up on a polyelectrolyte in a single dimension or spread over a two-dimensional surface.^12–14 However, Manning's theory, by virtue of relying exclusively on the comparison to the thermal energy level and a critical charge density, fails under conditions of elevated ionic strength or charge densities below the theorized onset of counterion condensation.^15,16 Several intricate accounts have detailed the influence of the size of a spherical particle and its curvature on the extent of the dissociation,^17–19 while concurrent studies focus on the reduced bulk ionic strength to renormalize the expected, bare surface charge to take counterion condensation into account.^20–22 The latter studies managed to provide methodologies to accurately describe the results of Monte-Carlo simulations of (theoretical) dissociating polyelectrolytes and nanoparticles and provide credible values for the decrease in surface charge in response to ionic strength independent of Manning's deliberations.

Yet the dissociation of individual, tightly packed charged moieties on a nanoparticle surface has never been scrutinized with regard to the chemical equilibrium of the individual dissociation reactions. The key consequence of the surface charge affecting a shift of the dissociation constants is that the acidity of the surface groups decreases significantly, which is reflected in potentiometric titrations.^22–24 Additionally, as the counterion condensation effect is not limited to protolysis reactions, other cations will condense onto the particle surface just as well.^4,5,13 This shift in the dissociation equilibria results in reduced ionic strength in the dispersion medium, which can be quantified by conductometry.^7,25

As such, both titrimetric analyses, despite being among the most common techniques to quantify surface charge, are tainted by fluctuating surface potentials, which, if not taken into account, may lead to misguided interpretations which are quantitatively significant.^26–30 Here, we suggest an alternative approach to determining the equivalence points in conductometric titration by quantifying the fraction of counterions entering the bulk solution to complement observations from conductivity or pH. The approach is based on monitoring immobilized, chargeable functional groups on a nanorod surface that can be quantitatively addressed by counterion condensation and Poisson–Boltzmann theories, leading to an updated dissociation equilibrium.

In order to prove the point in our case study, we have selected a challenging substrate to work with. Cellulose nanocrystals (CNCs) are essentially rigid, highly crystalline nanorods, but their top-down isolation from bio-based plant fibers means that slightly heterogenous morphologies with undefined cross-sections are encountered.^31,32 Nevertheless, they may be regarded as one-dimensional, roughly cylindrical particles with defined aspect ratios.^33,34 Our CNCs are stabilized in dispersion by phosphate moieties introduced by surface esterification, each of which carry one strongly acidic and one weakly acidic proton.³⁵ As such, their deprotonation may result in a potentially divalent charge, immobilized exclusively on the particle surface. The size distributions notwithstanding, we show how diligent application of existing physico-chemical theories^2,12,36–38 is sufficient to describe the dissociation phenomena accurately enough to simulate conductometric and potentiometric titration curves with good correlation and extract accurate values for surface charges and de facto dissociation constants.

Materials and methods

Materials

pCNC (1 mmol g⁻¹ phosphate as determined by elemental analysis; number average length of 116 nm and number average width and height of 6.25 nm; see, ESI†) were produced as outlined in our previous publication.³⁵ Concentrated hydrochloric acid solution (HCl, 34–37% assay, CAS 7647-01-0, Romil Chemicals Ltd), dilute hydrochloric acid solution (HCl, 0.1 M, CAS 7647-01-0, Titripur Reag., Merck), sodium hydroxide solution (NaOH, 1 M, CAS 1310-73-2, Titripur Reag., Merck) sodium chloride (NaCl, >99%, CAS 7740-23-5, VWR chemicals), sodium monophosphate dihydrate (NaH₂PO₄·2H₂O, CAS 13472-35-0, 98%, Supelco), ISO 17034-accredited AVS Titrinorm buffers (phthalate at pH 4.00, product number 32095-297; phosphate at pH 7.00, product number 32096-291; carbonate at pH 10.00, product number 32040-298, VWR chemicals) and Milli-Q water (18.2 MΩ cm resistivity) and conductivity calibrator solution (KCl in water, 1000 μS cm⁻¹ at 25 °C and 904 μS cm⁻¹ at 20 °C, product number 3167, YSI) were used without further purification.

Sodium hydroxide solution (1 M) was diluted to a concentration of 0.1 M with degassed Milli-Q water and used for conductometric titrations. Sodium chloride and for use as background electrolyte in conductometric titrations was dissolved in degassed Milli-Q water at 0.5 M.

pCNC synthesis

Phosphorylated cellulose nanocrystals (pCNC) were produced by the urea phosphate method as published previously³⁵ from cotton linters (Whatman 1 filter paper). The modified pulp was hydrolyzed with HCl gas immediately following the modification step, and then neutralized with NaOH (1 M).

We found that the subsequent washing procedure defines the obtained titration curves through the remaining counterions. Therefore, we include the treatment for our relevant batch of pCNC here. After HCl gas hydrolysis, wetting and neutralization (pH 7), the obtained modified and hydrolyzed pulp was washed exhaustively by centrifugation and redispersion in deionized water, until a pH of 8.9 and a conductivity of <20 μS cm⁻¹ of the supernatant were reached. The obtained pulp was then protonated by adjusting the pH to 1 using 37 wt% HCl solution. This solution at pH 1 was again treated to centrifugation and redispersion cycles until a final pH of 4.5 and a conductivity of <20 μS cm⁻¹ of the supernatant were reached. At this point, the individual nanocrystals were dispersed by microfluidization using a Microfluidics M-110P microfluidizer. The washed pulp was suspended in water at a concentration of 1 wt% and passed three times at 1500 bar through a pair of Z-type collision chambers with dimensions of 400 and 200 μm, respectively.

Unexpectedly, we found that this dispersion step has no noticeable consequence on the titration curves of 300 mg cellulose analyte in 500 mL water, although the surface area of the analyte was undoubtedly increased upon dispersion. Pulp, pCNC dispersed at pH 4.5 and pCNC dispersed at pH 7 yielded nearly identical datasets after normalization.

Conductometric titration

Conductometric titrations were carried out on the washed pulp using a 751 GDP Titrino dosing unit and a 712 Conductometer equipped with a conductivity measuring cell with Pt1000 thermocouple (Product number 6.0912.110) by Metrohm, following the protocol described by Ghanadpour et al.²⁵ Generally, equivalents of 300 mg dry, modified cellulose dispersed in water were added to 500 mL degassed Milli-Q water and 0.5 mL 0.5 M NaCl solution. The mixture was acidified with 5 mL 0.1 M HCl solution and titrated with 20 mL 0.1 M sodium hydroxide solution at 0.1 mL min⁻¹.

The electrode was calibrated using YSI conductivity standard solution.

Potentiometric titration

For more precise potentiometry measurements, where increased concentrations benefit the resolution of the measured curves, an equivalent of 860 mg modified cellulose was dispersed in 500 mL degassed Milli-Q water and acidified with 8.18 mL 0.1 M HCl solution. Correspondingly, comparative measurements were conducted on 8.60 mL 0.1 M sodium dihydrogen phosphate solution in 500 mL degassed Milli-Q water, acidified with 8.56 mL 0.1 M HCl solution. The nanocellulose dispersion was titrated with 40 mL 0.1 M NaOH solution at 0.1 mL min⁻¹, while the sodium dihydrogen phosphate solution, containing 1 eq. of sodium already, was titrated with 32 mL 0.1 M NaOH solution at 0.1mL min⁻¹. Simultaneously, pH-values were recorded for each increment using a Mettler Toledo SevenEasy pH-meter with an InLab Expert NTC30 electrode. We understand that the electrode is not designed for this prolonged use in low ionic strenght environments, but the reading obtained for the standard buffer solutions used for calibration (pH 4.00, 7.00 and 10.00) immediately after the titration revealed no discrepancies before and after the titration.

Titration curves of pulp (before) and dispersed pCNC (after microfluidization) were found to be identical. The dispersion state did not noticeably affect the neutralization process as observed by conductometry.

Impact of carbonation

As we conduct our experiments in ambient conditions, carbon dioxide from the air can dissipate into the analytes over time and convert to carbonic acid and carbonate over time, which impacts the neutralization reaction. As we will discuss in the next chapter, this is not entirely undesired, even though it creates a challenge for long titrations. In order to characterize the speed of this carbonation reaction, we conducted long term background measurements. In such experiments, 8 mL of NaOH titrant were dosed into degassed Milli-Q water at 0.1 mL min⁻¹, after which the conductivity of the solution was recorded overnight. By associating the steady conductivity decrease to the carbonation reaction progressing at constant speed, we estimated the impact of the carbonation reaction as described in the ESI.†

Normalization of the obtained datasets

Measured values depend on the applied voltage in the measurement as well as the electrode area and distance as illustrated in Fig. 2 and described in eqn (1), as well as the concentration of electrolytes in the volume between the electrode pair, which results in current flowing through the analyte. While the total amount of all electrolytes in the analyte is controlled, the volume may vary slightly from sample to sample due to varying dilutions and differing amounts of titrant volume. The measured conductivities are therefore specific to the individual analytes. In order to compare the separate experiments among each other, normalization is performed by multiplying the measured conductivities with the current volume of the analyte according to eqn (2).Here, σ is the measured conductivity, G is the conductance derived from current I and voltage U, d_el is the distance of the two electrode plates with area A_el, resulting in a volume of analyte between the two plates of V_el. V₀ is the initial volume of the analyte and V_B is the volume of added titrant.

Fig. 2 Schematic representation of the measurement geometry of the electrode. The current I transported between the two electrodes with area A_el and distance d_el at an applied voltage U yield conductance G and conductivity σ.

The physical meaning of the conversion from conductance to conductivity in eqn (1) is to scale the measured raw conductance by the measurement geometry. Our normalization according to eqn (2) extends this frame of the raw measurement to the entirety of the analyte volume. The normalized values are in between conductivity and conductance and can be converted back to conductivities if divided by a reference volume of the analyte or give conductances if divided by the (constant) square of the electrode distance. Nevertheless, in the following, we will refer to this normalized quantity as conductance.

Secondly, in order to compare the individual datasets containing slightly varying amounts of analyte and excess HCl, the x-axis was normalized to the amount of phosphate surface groups as determined by elemental analysis for the cellulosic analytes; correspondingly the sodium phosphate salt titration curve was normalized to the amount of phosphate in the analyte. This step allows us to account for the slight contamination through the carbonation reaction as well. The volumes of added NaOH were converted into equivalents according to eqn (3) for phosphorylated nanocellulose and eqn (4) for sodium phosphate standards.

Here, eq. refers to the number of equivalents of base added to the system, n_base is the amount of NaOH added in the titration, n_HCL is the amount of HCl added to initially acidify the analyte, n_Phosphate is the amount of phosphate in the analyte as determined by elemental analysis, and n_carbonate is the amount of carbonate that entered the analyte at the considered point in the experiment. V_B, V_HCL and V_NAH2PO4 refer to the volumina of base, HCl solution and sodium dihydrogen phosphate standard solutions of the respective concentrations c_B, c_HCL, and c_NAH2PO4. Furthermore, and are the perceived concentrations of the NaOH solution and the air affecting carbonation as outlined in the ESI.†

We emphasize that m_analyte and c_P,analyte refer to the dry mass of cellulose in the aliquot of cellulose dispersion added to the analyte solution and the concentration of phosphate groups in the (dry) cellulosic material. We also consider the stoichiometric amounts of carbonate rather than its valency as the neutralization of bicarbonate occurs at elevated pH and does not influence the first two equivalence points of the phosphate neutralization that we are interested in. The surface area of the electrolyte was assumed to be constant throughout the measurements.

Results and discussion

We conducted conductometric titrations to assess the diffusion of the ions in the bulk phase, which are mainly due to protons, sodium, chloride and hydroxide ions, and potentiometric titrations, which detect the bulk concentrations of protons and hydroxide ions exclusively. Potentiometric and conductometric titration curves of phosphorylated cellulose nanocrystals and an equivalent amount of sodium phosphate are shown in Fig. 3.

Fig. 3 Conductometric and potentiometric titration curves for phosphorylated CNC (black and red), respectively and an equivalent amount of monosodium phosphate (grey and pink, respectively). The curves for monosodium phosphate are shifted along the x-axis to reflect the initial partial neutralization.

Despite containing the same amounts of phosphate units, the conductivity of the CNC-containing analyte is significantly reduced. In part, this is due to the discrepancies in molar conductivities of the ionic and immobilized phosphates. The major cause for the discrepancy, however, is the impact of counterion condensation,⁷ which significantly reduces the ionic strength of the bulk phase around the colloidal analyte. Simultaneously, the comparison of the pH curves reveals the perceived increase in basicity of the bound phosphate moieties in accordance with previously reported studies.^18,24 Also, the individual shapes of the conductometric titration curves of immobilized^27,28,39 and ionic²⁵ phosphates have been reported previously, but the peculiar shape of the former ones went without more detailed discussion. It is evident in Fig. 3 and demonstrated in Fig. S3 (ESI†), however, that the analysis of the immobilized phosphate content by linearization leads to a significant overestimation of the actual result. In order to eventually address this issue, it is necessary to understand the nature and extent of the individual ionic dissociation processes that are the foundation of the analytic method.

The dissociation reactions connecting the different phosphate surface group species are shown in Fig. 4. We are neglecting the formation of disodium and trisodium orthophosphate here, given that, for one, no equilibrium constants to describe these processes accurately were available and, for another, their concentrations can be considered negligible in high dilution and below pH 12. Besides the phosphate analyte, pH-active carbonate species are present in the system. These arise from the limited exposure of the analyte solution to ambient air, allowing for carbonation reactions to take place. The species we deem necessary to be considered and their dissociation processes are shown in Fig. 5. Again, while it is theoretically possible that disodium carbonate is formed in solution, under the highly diluted conditions we deem its concentrations negligible.

Fig. 4 AFM image of the analyzed pCNC (a) and concurrent dissociation reactions of phosphate moieties involving protons and sodium ions (b). R refers to the cellulose entity the surface groups are anchored to or, in the case of free orthophosphoric acid, R refers to the least acidic proton.

Fig. 5 Carbonation reaction responsible for the formation of bicarbonate and carbonate, which consume titrant and contribute to the conductivity of the analyte.

All individual reactions are described by their own individual equilibrium constants, which are summarized in Table 1. The equilibrium constant for the interconversion of sodium phosphate and disodium phosphate was estimated as shown in the ESI.†

Table 1 Equilibrium constants for the dissociation or association processes j of the relevant phosphate and carbonate species. Bold values are adapted from the indicated reference, the other values were computed on that basis

j	Dissociation constant		Association constant		Ref.
	Kdis	pKdis	Kas	pKas
a1	1.10 × 10^−2	1.96	91.2	−1.96	40
a2	6.17 × 10^−8	7.21	1.62 × 10^−7	−7.21	40
a3, eff	6.31 × 10^−13	12.2	1.28 × 10^−12	−12.2	40
d1	0.56	0.25	1.78	−0.25	41
d2	2.09 × 10^−12	1.68	47.9	−1.68	41
d3	0.22	0.65	4.47	−0.65	42 Estimate (see ESI)
a4, eff	4.47 × 10^−7	6.36	2.24 × 10^−6	−6.36	43
a5	4.68 × 10^−11	10.35	2.14 × 10^−10	−10.35	43
d4	1.51	−0.18	0.66	0.18	43
d5	9.77 × 10^−2	1.01	10.2	−1.01	43
Kw	1.0 × 10^−14	14	1.0 × 10^−14	−14

---

With the knowledge of the accurate volume of the considered analyte and the total amounts of all components in the analyte, the concentrations of all individual species in the interconnected dissociation equilibria can be calculated. The details will be outlined in the following.

Firstly, the total amounts of chloride as well as sodium, carbonates, and phosphates are defined by the mass balance of all relevant species. Trivially, the concentration of chloride ions is equal to the total concentration of chloride in the system (eqn (5)). Analogously, however, the total amount of all species containing sodium (eqn (6)), phosphate (eqn (7)) and carbonate (eqn (8)) is known and can be expressed as follows:

[Σ(chloride)] = [Cl⁻] (5)

[Σ(sodium)] = [Na⁺] + [NaRPO₄⁻] + [HNaRPO₄] + 2[Na₂RPO₄] + [NaHCO₃] + [NaCO₃⁻] (6)

[Σ(phosphate)] = [H₂RPO₄] + [HRPO₄⁻] + [RPO₄²⁻] + [NaHRPO₄] + [NaRPO₄⁻] + [Na₂RPO₄] (7)

[Σ(carbonate)] = [CO_2,aq] + [HCO₃⁻] + [NaHCO₃] + [NaCO₃⁻] + [CO₃²⁻] (8)

Here, [i] refers to the concentration of compound i, while [Σ] refers to the sum of the individual concentrations of interrelated ions.

The exact concentrations of all individual species contained in eqn (6)–(8) can be assessed when all dissociation processes are at their equilibrium states. At equilibrium, the ratios of dissociated and associated species are defined by the individual dissociation constants, provided in Table 1. However, as discussed in the introduction, in the presence of surface potentials, the equilibrium constants for the dissociation of the phosphate surface groups require additional correction. If and when this phenomenon becomes relevant can be assessed by considering the electrostatic potential energy between two charges as expressed in eqn (9)

with

In eqn (9), (10) and (11), U is the electrostatic potential energy, ε₀ and ε_r are the vacuum permittivity and the relative permittivity of water, q₁ and q₂ are the amounts of charge interacting, and r₁ and r₂ are the locations of two interacting charges, with d_charge being their Euclidian distance, k_B is the Boltzmann constant, T is the absolute temperature and λ_B is the Bjerrum length. If the distance between the two charges is large, their electrostatic interaction is negligible. In close proximity, however, the electrostatic potential energy quickly increases as the individual electric fields overlap. According to Manning, the threshold for counterion condensation is reached, when the electrostatic potential between two neighboring charges non longer decreases below the thermal energy level according to eqn (10). In this context, it is sensible to express the distance between two charges in terms of the Bjerrum length λ_B (see, eqn (11)), which defines the distance, in which the potential of a single elementary charge decreases to thermal level. Following Manning's theory, if the charges are separated by less than a Bjerrum length, the superposition of their individual potentials leads to significantly enhanced electrostatic interactions that induce counterion condensation. This is an absolute criterion in Manning's original theory, which does not allow for degrees of surface interactions. Either the threshold is exceeded or not, which makes this theory unsuited for low surface charges. Similarly, by taking into account chemical potentials of the dissociating ions, Manning was able to describe selectivities for mixtures of counterions,⁴⁴ but his model fails, for instance, for nearly quantitative protonation of polyacids.¹⁶ Nevertheless, the concept of separately considering the effects of the dissociation free energies and the electrostatic contributions on the overall dissociation reactions, as outlined in Manning's later work,¹² can be applied more universally, to overcome the aforementioned limitations. This is illustrated in Fig. 1d and can be expressed according to eqn (12) and (13).

ΔG⁰_total = ΔG_{dissociation0} + U_site (12)

By treating the electrostatic interactions of the site with the surrounding surface as an additive component to the free energy ΔG of the main dissociation reaction, it follows that this contribution can be applied to the equilibrium constant as a correctional factor K_condensation according to eqn (13). This simplification means that we neglect entropic differences between the colloidal and ionic phosphates that originate from the presence of an electrostatic double layer around the colloidal particles, as well as potential changes in the hydration of immobilized and dissociated phosphate moieties. However, there is precedent for the separate contemplation of the ionic interactions between the two dissociating charges and a background potential,^2,5,45 and treating entropic changes as negligible seems reasonable when considering that for both ionic and immobilized phosphate, dissociation events increase the total number of particles by the same increment: 1. this means that the absolute entropies of ionic and immobilized phosphate systems must be differing due to the anionic and cationic charge concentrations around the particles, but the entropy changes due to individual dissociation events can be considered nearly identical.

We have detailed our assessment of the electrostatic potential energy and its impact on the dissociation equilibria in Appendix A (ESI†). These deliberations are equivalent to previously published charge renormalization methodologies.^20,21,46,47 We, too, consider the promotion of ion association to the surface, which reduces the degree of dissociation of the potentially charged sites, but from the perspective of the individual dissociation events. This means that our terminology varies from that of established charge renormalization considerations. Charge normalization is founded in the solution of the nonlinear Poisson–Boltzmann equation describing the charge distribution and potential decrease within the confines of a Wigner–Seitz cell around a colloidal particle with a defined bare charge. With the knowledge of the electrostatic potential at the cell boundary and the bulk ionic strength, it is possible to linearize the Poisson–Boltzmann equation and calculate effective surface charges, by considering the integrated charge density and balancing it against the net charge of the Wigner–Seitz cell.^20,46 In this context, the bare surface charge is considered partially neutralized by unspecific ion association, which neglects ionic affinities. Nevertheless, accurate results are obtained when the association constants of competing counterions are several magnitudes apart.^22,48 Here, instead, we consider the surface dissociation processes as directly impacted by the surface potential, meaning phosphate functionalities carrying bound protons or forming ion pairs with sodium ions are considered as uncharged. Therefore, in our consideration, we follow the affinities of individual counterions for individual surface functionalities through their corresponding equilibrium constants, and emerging surface charges are directly comparable to the otherwise renormalized values.²² This methodology allows us to track both strong and weak acidic functionalities, as well as associates of both protons and sodium ions, simultaneously.

Following eqn (13), the influence of the surface energy can be expressed as an additional constant K_condensation. As such, eqn (14) and (15) describe the deprotonation of phosphate surface groups, while eqn (16)–(18) describe the analogous dissociation of sodium phosphate entities and sodium hydrogen phosphate.

Here, K_j refers to the ionic dissociation constant of reaction j, using the same terminology as in Fig. 4 and 5 as well as Table 1.

Analogous expressions describe the protonation (eqn (19) and (20)) and the sodium dissociation equilibria (eqn (21) and (22)) for bicarbonate and carbonate, respectively. As these are dissolved, ionic species, counterion condensation plays no role here.

These equilibrium conditions yield the ratios of the considered species if the concentration of free protons is known. This concentration, in turn, follows from the charge balance of all species in the system (eqn (23)) and the ion product of water (eqn (24)).

0 = ([OH⁻] + [Cl⁻] + [HRPO₄⁻] + 2[RPO₄²⁻] + [NaRPO₄⁻] + [HCO₃⁻] + 2[CO₃²⁻] + [NaCO₃⁻]) − ([H⁻] + [Na⁺]) (23)

By rearranging eqn (5)–(24) into matrix form as shown in Appendix B (ESI†), we obtain a system of equations described by eqn (25), that defines the equilibrium concentrations of all dissociation equilibria considered here. In this context, it is convenient to convert the expressions involving dissociation constants into logarithmic form. This means that the concentration quotients are converted into sums, which are easier to handle, and which prevent the consideration of negative values for the real concentrations. The logarithmic values of the concentrations can be obtained by numerically solving eqn (25).

C_ln(V_B) = M⁻¹ × K(V_B) (25)

Here, C_ln(V_B) is the vector of logarithmic concentrations at titrant volume V_B, as defined in eqn (B13) and (B17) (ESI†) for immobilized and ionic phosphate, respectively, M⁻¹ is the inverse of the matrix M as defined in eqn (B11) and (B15) (ESI†) for immobilized and ionic phosphate, respectively. K(V_B) is the vector containing the system information at titrant volume V_B as defined in eqn (B12) and (B16) (ESI†) for immobilized and ionic phosphate, respectively.

The concentrations obtained by eqn (25) describe the dissociation of the surface groups accurately if the estimate account for the counterion condensation is accurate. However, the bulk concentrations of the mobile ions, which account for the measured conductances, remain veiled by the diffuse layer of counterions around the surface charge. As illustrated in Fig. 6 and 7, the electric surface potential attracts counterions around the particle, compensating the charge.^37,49 Fig. 6 also emphasizes the fact that the zeta potential – as the potential along the slipping plane – depends on several external factors and is not directly proportional to the surface charge density. For reference, we have reported the zeta potential values of pCNCs in our earlier publication.³⁵

Fig. 6 Schematic representation of the electrostatic potential decrease throughout the diffuse layer, situated between bulk phase (r > r_B) and Stern layer (r_P < r < r_P + r_ion). The measured zeta potentials (yellow) tend to be slightly larger than thermal potential,³⁵ placing the slipping plane closer to the particle surface than our assumed boundary plane.

Fig. 7 Schematic representation of the distribution of cations (red) and anions (blue) around the cross section of a cylindrical cellulose particle with negative surface charge. The dashed lines indicate the separation between Stern layer and Gouy–Chapman layer (red), the slipping plane (yellow) and the boundary plane (blue), which separates the Gouy–Chapman layer from the bulk phase.

The loosely bound ions close to the surface are kinetically tied to the particles, which reduces their bulk activities and their ability to conduct electricity. Using the known surface charge, it is possible to evaluate the concentrations of cations and anions in the diffuse layer by applying Poisson–Boltzmann statistics,^36,50,51 as has been demonstrated by studies predicting the osmotic pressure of colloidal dispersions.^20,21

By integrating the ion densities predicted by the Poisson Boltzmann approach over the volume of the electrostatic double layer, it is possible to “count” the ions that are immobilized close to the particle surface.^20,36 Thus, following the methodology outlined in Appendix C (ESI†), we can obtain activity values for the dissociated ions calculated in eqn (25) by balancing the number of ions immobilized in the diffuse layer against the total number of ions in the system. These activity values can then be applied to obtain the bulk concentrations of the individual species according to eqn (26) and (27).

c^bulk(V_B) = γ(V_B)·c(V_B) (26)

C^bulk_ln(V_B) = C_ln(V_B) + γ_ln(V_B) (27)

Here, c(V_B) is the vector containing the real overall concentrations of dissociated ions derived from C_ln(V_B) by raising e to the individual vector elements’ power. Correspondingly, bulk concentrations c^bulk(V_B) are derived by element-wise multiplication of the overall concentrations of dissociated ions with the bulk activities γ(V_B). The logarithmic equivalent based on C^bulk_ln(V_B), C_ln(V_B) and γ_ln(V_B) is simply an addition, which is easier to compute.

The real (non-logarithmic) bulk concentrations, in turn, can be converted into conductivities according to eqn (28) using the individual limiting equivalent conductivities Λ⁰_i listed in Table 2.

Table 2 Equivalent conductivities of the relevant ionic species bold values are adapted from literature, the conductivity of the particles is estimated in Appendix D (ESI) and the remaining conductivities are estimated based on 2% °C⁻¹ estimates

Λ^0i	Literature values at 25 °C	Values at 20 °C	Ref.
a In literature, the conductivities of tribasic monosodium phosphate and monosodium carbonate ions have been estimated to be equal to hydrogen phosphate and bicarbonate ions, respectively.^41,54
	349.85 (mS cm^−1) (L mol^−1)	325.52 (mS cm^−1) (L mol^−1)	52
	198.3 (mS cm^−1) (L mol^−1)	179.2 (mS cm^−1) (L mol^−1)	52
	50.1 (mS cm^−1) (L mol^−1)	44.8 (mS cm^−1) (L mol^−1)	52
	76.4 (mS cm^−1) (L mol^−1)	69.7 (mS cm^−1) (L mol^−1)	52
	32 (mS cm^−1) (L mol^−1)	28.43 (mS cm^−1) (L mol^−1)	52
	57.2 (mS cm^−1) (L mol^−1)	50.82 (mS cm^−1) (L mol^−1)	52
	^a	^a	41
	46 (mS cm^−1) (L mol^−1)	41.4 (mS cm^−1) (L mol^−1)	53
	^a	^a	54
	69.3 (mS cm^−1) (L mol^−1)	62.4 (mS cm^−1) (L mol^−1)	55
Λ^0CNC	Eqn (D8) (ESI)

---

The conductivity of the CNC depends on their charge as well as the mobility of the counterions in the double layer and can be estimated as described in Appendix D (ESI†).

The conductivity data that can be obtained with these calculations are compared to the experimental data in Fig. 8a, while the measured data for the pH during the titration are compared to the calculated values in Fig. 8b.

Fig. 8 Simulated conductivity (A) and potentiometry (B) curves for phosphate surface groups (blue) and dissolved phosphate (red) compared to the measured datasets (black) and calculated logarithmic values for the condensation constant (pK_condensation) impeding the dissociation (dashed blue line).

It is clear from Fig. 8a and b that for dissolved sodium dihydrogen phosphate, the calculated data for both pH and conductivity are in excellent agreement with the measured data. Evidently, our model for the dissociation equilibria of the phosphate species describes the observed phenomena accurately. By contrast, both calculated datasets for the immobilized phosphate surface groups accurately describe the trends in the inhibition of the dissociation of the surface groups which leads to higher observed pK_a values as well as lower overall conductance and widened curve shapes in conductometry, but the clear discrepancies in theoretical and practical values expose flaws in our model. As such, the conductance curve in Fig. 8 is generally underestimating the conductance until sodium hydroxide is added in excess, raising the bulk anion concentration in particular and ionic strenght in general. This suggests that the corrections we employed for the counterion condensation effect are too strong. Specifically, this means that either the surface potential or the amount of ions trapped in the electrostatic double layer are overestimated. Additional experiments in environments with different background ionic strenghts, which are shown in Fig. 9, indicate, that it is in fact the overestimation of the surface potential that is to blame.

Fig. 9 Conductometric titration curves of immobilized phosphate in the presence of varying amounts of background electrolyte. The analyte dispersions contain 0.4 mM (0.66 eq., black), 1.0 mM (1.66 eq., red), 1.6 mM (2.66 eq., green), and 2.2mM (3.66 eq., blue) Chloride as inert background electrolyte, originating from HCl and NaCl.

The curves shown in Fig. 9 show similar discrepancies between modeled and measured datasets. The varying ionic strengths have a significant impact on the Debye length, which describes the screening of the surface charges as well as the thickness of the electrostatic double layer. Despite the increased screening at high ionic strengths, however, the observed discrepancies remain similar.

Consulting the pH curves in Fig. 8b further supports these deliberations. It is evident, that the calculated pH values are consistently larger than the observed values in the acidic regime, illustrating the overly reduced acidity in the theoretical dataset. In the basic regime, in turn, the pH is underestimated, which can be attributed to the growing concentration of sodium ions in the system. Our model assumes that too many sodium ions are associated to the surface groups, facilitating the dissociation of protons. It must be considered, however, that past the first equivalence point, the discrepancies between measured and calculated pH describe discrepancies in hydroxide concentration that are lower than 0.01% and well below the uncertainties of the concentrations of the standard solutions used in the titrations.

Possible causes for the overestimation of the surface potentials lie in the total amounts of surface area and charge, respectively, as well as in their arrangement on the surface and the strength of their interaction. While we can regard the overall charge content as reliable, having confirmed it via elemental analysis and found corresponding results from the ash contents in thermogravimetric analyses, the size of the crystallites is less clear. The CNC have been isolated in a top-down approach from natural fibers and therefore display distributions around median values for length and, crucially, height,^33,35,56 which translates to the surface area of the individual particles. Furthermore, while it is clear from their dimensions, that the particles are aggregates of several crystallites, the actual arrangement of these elemental fibril fragments^57,58 and that of cellulose chains along their cross-sections⁵⁹ remains unclear. Therefore, while we cannot specify where exactly the surface groups are located on the particle surface, the most prudent assumption is that of cylindrical particles of average diameter with a homogenous charge distribution on the mantle area as described in Appendix A (ESI†).

Despite this simplification, it is evident that our model is still a good approximation of the observed data. It correctly describes the significant decrease in conductivity compared to the dissolved phosphate, as well as the emergence of a local minimum close to the first equivalence point ahead of the plateau phase leading up to the second equivalence point. Around the second equivalence point, our model predicts an idealized linear behavior that is not reflected in the measured data, but the final increase in conductance due to excess sodium hydroxide is reflected very accurately. Therefore, while reviewing them critically, the calculated data can still give insight into the processes that need to be considered in the analysis. First and foremost, the counterion condensation parameter, which is displayed in logarithmic form in Fig. 8b, illustrates that the acidic strength of the phosphate surface groups is decreased significantly, raising the observed pK_a-value for the first acidic dissociation reaction from 3.5 to 4.5 depending on the surface charge. Consequently, during the initial deprotonation, the surface groups are strongly acidic as are dissolved phosphate ions, but then the acidic strength of the surface groups decreases significantly. That is the cause for both the significant broadening of the titration curve and the decreased initial conductivity. Simultaneously, the model estimate of effective pK_a values in the order of 3.5 to 4.5 for the stronger acidic proton illustrates that the analysis of the conductometric titration curve merely by linearization does not apply here. The surface groups are transitioning from strongly acidic to weakly acidic as a surface potential is generated, which precludes accurate quantification.

Instead, the minimum in the curves preceding the plateau phase is a more reliable feature for identifying the equivalence point. This claim is supported by the change in the bulk sodium concentration according to our model, which is shown in Fig. 10.

Fig. 10 Fraction of sodium ions entering the bulk solution with regard to added titrant in the case of immobilized phosphate (blue) and dissolved phosphate (red) (a), and illustration of the capture of sodium ions by the immobilized phosphate groups as opposed to ionic phosphate (b).

It is clear from Fig. 10, that in the case of dissolved phosphate ions, the association of sodium ions is negligible before the first equivalence point is reached. Subsequently, the formation of divalent hydrogen phosphate ions leads to the simultaneous formation of sodium hydrogen phosphate associates, which extends past the second equivalence point, when the sodium ion concentration in the system keeps increasing. Therefore, the amount of sodium ions in the bulk phase increases, but a fraction of the added ions associates to phosphate moieties. This effect is substantially more significant for the immobilized phosphate surface groups. In fact, the amount of sodium ions entering the bulk phase of the solution between the first and second equivalence points is near zero, illustrating that a maximal surface charge is reached. At this point, further deprotonation is accompanied with immediate association of a sodium ion. However, at the first equivalence point, the rate of sodium entering the bulk phase of the system spikes. This is an indirect consequence of the decreased acidity of the surface groups, shifting the jump in pH at the first equivalence point from 3–5 to 4–6 (see, Fig. 8b). This means that in the case of the immobilized phosphate units, any carbonic acid impurities in the system are being neutralized inbetween the first and second protons of the phosphate units, immediately after the first equivalence point is reached. During the deprotonation of the carbonic acid, the surface charge of the cellulose fibers increases to a lower degree, meaning that more free sodium ions enter the bulk phase. This increase in carbonate and sodium ions in the bulk phase is visible in the titration curves in Fig. 8a and Fig. 9 as the slight increase in conductivity after the first equivalence point. By contrast, in the case of dissolved phosphate, where the deprotonation of the dihydrogen phosphate ions starts already at a pH level of 5, both dihydrogen phosphate (pK_a 7.21) and carbonic acid (apparent pK_a 6.35, see, Table 1) deprotonate simultaneously when the pH increases.

Lastly, our modeled data give access to surface charge values, which take into account the incomplete dissociation and the screening of the surface charges by the surrounding diffuse layer. We should point out here, that these surface charge values correspond to the renormalized values discussed in other studies.^{20–22,46,47} The net charges at the boundary layer are not related to those studies, as they do not consider this boundary between diffuse layer and bulk phase. Our data are shown in Fig. 11.

Fig. 11 Calculated values for the specific surface charge of the CNC (left) and net surface charge after screening by the diffuse layer (right). The values are calculated for analyte dispersions containing 0.4 mM (0.66 eq., black), 1.0 mM (1.66 eq., red), 1.6 mM (2.66 eq., green), and 2.2 mM (3.66 eq., blue) chloride background electrolyte, the yellow curve is calculated for the titration of a triple amount of phosphate analyte with 1.66 mM (0.95 eq.) chloride.

As shown in Fig. 8 and 9, the calculated values are in reasonable agreement with the measured data and accurately describe the decrease in conductivity for the particles caused by incomplete dissociation. This means that the calculated surface charge densities of less than 0.5 mmol g⁻¹, which translates to less than 0.35 mmol g⁻¹ at the boundary plane need to be taken seriously, even if this suggests that less than 25% of potentially chargeable sites on the surface are actually dissociated at any given time. It must be noted that our model describes a decrease in surface charge leading up to the second equivalence point. This is where calculated and measured data have their largest discrepancy (see, Fig. 8). Therefore, this reduction in surface charge is most likely an artifact caused by the overestimation of the counterion condensation effect. Furthermore, our model predicts an increase in surface charge with increasing ionic strength. While this seems counterintuitive, given more efficient screening of the surface charge in conditions with reduced Debye lengths, the cause for this effect can be seen in eqn (A4) (ESI†). The surface potential also depends on the Debye length, meaning that an increased charge density in the diffuse layer also causes screening of the surface charges’ interactions among each other. This means that the effective surface potential is reduced with increasing ionic strength, shifting the surface equilibrium reactions back towards dissociation. The dependance of the surface charge on the ionic strength is further illustrated in Fig. 12.

Fig. 12 Dependance of the surface charge on the ionic strength.

During the neutralization processes, the bulk ionic strength remains nearly constant, the net reaction being the replacement of protons in the system with sodium ions. Only after the second equivalence point is reached, the ionic strength increases due to excess sodium and hydroxide ions accumulating in the system. In this context, Fig. 12 illustrates that our model calculates very similar surface charges for comparable ionic strengths for all experiments, regardless of the actual particle concentrations.

Conclusions

We have demonstrated the significance of counterion condensation for the dissociation of phosphorylated cellulose nanocrystals. This effect is the root cause of the non-ideal shape of both conductometric and potentiometric titration curves, which have been reported before, but, to our knowledge, never been discussed in detail. Here, we managed to develop mathematical descriptions for the surface potential, counterion condensation corrected dissociation constants and screening of the particles, enabling us to simulate the titration curves with reasonable success. Furthermore, we proposed an alternative method to deduce the extent of dissociation from conductometric data based on calculating the fraction of counterions entering the bulk solution instead of merely observing conductivity or pH. This methodology further enabled us to speciate the associated ions, which represents a refinement of previous works on charge renormalization. While not perfect, our model delivers insight into the specifics of the dissociation processes, allowing us to estimate the counterion condensation corrected surface charge at 0.2–0.5 mmol g⁻¹ for CNCs carrying 1 mmol g⁻¹ phosphate surface groups.

Author contributions

Marcel Kröger: conceptualization, methodology, formal analysis, visualization, investigation, writing – original draft, writing – review & editing. Kaarlo Nieminen: software, formal analysis, writing – review & editing. Han Tao: visualization, investigation, writing – review & editing. Eero Kontturi: conceptualization, writing – original draft, writing – review & editing, funding acquisition, project administration.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the ESI.†

Acknowledgements

We acknowledge FinnCERES Flagship funded by Research Council of Finland (grant no. 318890 and 318891).

References

1. J. Lyklema, H. P. van Leeuwen and M. Minor, Adv. Colloid Interface Sci., 1999, 83, 33–69.
2. P. H. Debye, Ernst, Phys. Z., 1923, 24, 185–206.
3. I. V. Pobelov, in Encyclopedia of Interfacial Chemistry, ed. K. Wandelt, Elsevier, Oxford, 2018, pp. 299–315.
4. J. R. Huizenga, P. F. Grieger and F. T. Wall, J. Am. Chem. Soc., 1950, 72, 4228–4232.
5. R. O. James, J. A. Davis and J. O. Leckie, J. Colloid Interface Sci., 1978, 65, 331–344.
6. M. E. Labib and A. A. Robertson, J. Colloid Interface Sci., 1980, 77, 151–161.
7. H. P. v Leeuwen, R. F. M. J. Cleven and P. Valenta, Pure Appl. Chem., 1991, 63, 1251–1268.
8. G. S. Manning, J. Phys. Chem., 1975, 79, 262–265.
9. G. S. Manning, J. Chem. Phys., 1969, 51, 924–933.
10. G. S. Manning, J. Chem. Phys., 1969, 51, 934–938.
11. G. S. Manning, J. Chem. Phys., 1969, 51, 3249–3252.
12. G. S. Manning, J. Phys. Chem. B, 2007, 111, 8554–8559.
13. J. A. Davis, R. O. James and J. O. Leckie, J. Colloid Interface Sci., 1978, 63, 480–499.
14. M. L. Henle, C. D. Santangelo, D. M. Patel and P. A. Pincus, Europhys. Lett., 2004, 66, 284.
15. E. Trizac and G. Téllez, Phys. Rev. Lett., 2006, 96, 038302.
16. L. Yan, Rep. Prog. Phys., 2002, 65, 1577.
17. H.-K. Tsao, Y.-J. Sheng and C.-Y. D. Lu, J. Chem. Phys., 2000, 113, 10304–10312.
18. D. Wang, R. J. Nap, I. Lagzi, B. Kowalczyk, S. Han, B. A. Grzybowski and I. Szleifer, J. Am. Chem. Soc., 2011, 133, 2192–2197.
19. D. A. J. Gillespie, J. E. Hallett, O. Elujoba, A. F. Che Hamzah, R. M. Richardson and P. Bartlett, Soft Matter, 2014, 10, 566–577.
20. E. Trizac, L. Bocquet, M. Aubouy and H. H. von Grünberg, Langmuir, 2003, 19, 4027–4033.
21. S. Pianegonda, E. Trizac and Y. Levin, J. Chem. Phys., 2007, 126.
22. A. Bakhshandeh, D. Frydel and Y. Levin, Langmuir, 2022, 38, 13963–13971.
23. J. Lützenkirchen, T. Preočanin, D. Kovačević, V. Tomišić, L. Lövgren and N. Kallay, Croat. Chem. Acta, 2012, 85, 391–417.
24. G. Charron, D. Hühn, A. Perrier, L. Cordier, C. J. Pickett, T. Nann and W. J. Parak, Langmuir, 2012, 28, 15141–15149.
25. I. M. Kolthoff, Ind. Eng. Chem. Anal. Ed., 1930, 2, 225–230.
26. M. Ghanadpour, F. Carosio, P. T. Larsson and L. Wågberg, Phosphorylated Cellulose Nanofibrils: A Renewable Nanomaterial for the Preparation of Intrinsically Flame-Retardant Materials, 2015, vol. 16, pp. 3399–3410 DOI:10.1021/acs.biomac.5b01117.
27. A. Ait Benhamou, Z. Kassab, M. Nadifiyine, M. H. Salim, H. Sehaqui, A. Moubarik and M. El Achaby, Cellulose, 2021, 28, 4625–4642.
28. Y. Noguchi, I. Homma and Y. Matsubara, Cellulose, 2017, 24, 1295–1305.
29. E.-H. Ablouh, Z. Kassab, F.-Z. Semlali Aouragh Hassani, M. El Achaby and H. Sehaqui, RSC Adv., 2022, 12, 1084–1094.
30. M. Zhao, S. Fujisawa and T. Saito, Distribution and Quantification of Diverse Functional Groups on Phosphorylated Nanocellulose Surfaces, 2021, vol. 22, pp. 5214–5222 DOI:10.1021/acs.biomac.1c01143.
31. T. Maloney, J. Phiri, A. Zitting, A. Paajanen, P. Penttilä and S. Ceccherini, Carbohydr. Polym., 2023, 319, 121166.
32. T. Benselfelt, N. Kummer, M. Nordenström, A. B. Fall, G. Nyström and L. Wågberg, ChemSusChem, 2023, 16, e202201955.
33. X. M. Dong, J.-F. Revol and D. G. Gray, Cellulose, 1998, 5, 19–32.
34. R. J. Moon, A. Martini, J. Nairn, J. Simonsen and J. Youngblood, Chem. Soc. Rev., 2011, 40, 3941–3994.
35. M. Kröger, O. Badara, T. Pääkkönen, I. Schlapp-Hackl, S. Hietala and E. Kontturi, Biomacromolecules, 2023, 24, 1318–1328.
36. D. Stigter, J. Colloid Interface Sci., 1975, 53, 296–306.
37. D. C. Grahame, Chem. Rev., 1947, 41, 441–501.
38. W. H. Lee and R. J. Wheaton, J. Chem. Soc., Faraday Trans. 2, 1979, 75, 1128–1145.
39. E.-H. Ablouh, F. Brouillette, M. Taourirte, H. Sehaqui, M. El Achaby and A. Belfkira, RSC Adv., 2021, 11, 24206–24216.
40. K. Schrödter, G. Bettermann, T. Staffel, F. Wahl, T. Klein and T. Hofmann, Ullmann's Encyclopedia of Industrial Chemistry, 2008 DOI:10.1002/14356007.a19_465.pub3.
41. A. D. Pethybridge, J. D. R. Talbot and W. A. House, J. Solution Chem., 2006, 35, 381–393.
42. P. G. Daniele, A. De Robertis, C. De Stefano, A. Gianguzza and S. Sammartano, J. Solution Chem., 1991, 20, 495–515.
43. A. Stefánsson, P. Bénézeth and J. Schott, Geochim. Cosmochim. Acta, 2013, 120, 600–611.
44. G. S. Manning, J. Phys. Chem., 1984, 88, 6654–6661.
45. O. Ciftja, A. Babineaux and N. Hafeez, Eur. J. Phys., 2009, 30, 623.
46. S. Alexander, P. M. Chaikin, P. Grant, G. J. Morales, P. Pincus and D. Hone, J. Chem. Phys., 1984, 80, 5776–5781.
47. E. Trizac and Y. Levin, Phys. Rev. E, 2004, 69, 031403.
48. G. H. Zimmerman and R. H. Wood, J. Solution Chem., 2002, 31, 995–1017.
49. O. Stern, J. Electrochem. Appl. Phys. Chem., 1924, 30, 508–516.
50. H. Saboorian-Jooybari and Z. Chen, Results Phys., 2019, 15, 102501.
51. J. A. Bertolotto and J. P. Umazano, J. Phys.: Condens. Matter, 2022, 34.
52. E. N. Tsurko, R. Neueder, J. Barthel and A. Apelblat, J. Solution Chem., 1999, 28, 973–999.
53. J. Kendall, J. Am. Chem. Soc., 1916, 38, 1480–1497.
54. O. Bernard and J. Aupiais, J. Mol. Liq., 2018, 272, 631–637.
55. C. B. Monk, J. Chem. Soc., 1949, 429–431.
56. T. Pääkkönen, P. Spiliopoulos, A. Knuts, K. Nieminen, L.-S. Johansson, E. Enqvist and E. Kontturi, React. Chem. Eng., 2018, 3, 312–318.
57. M. Martínez-Sanz, F. Pettolino, B. Flanagan, M. J. Gidley and E. P. Gilbert, Carbohydr. Polym., 2017, 175, 450–463.
58. K. Daicho, K. Kobayashi, S. Fujisawa and T. Saito, Angew. Chem., Int. Ed., 2021, 60, 24630–24636.
59. J. D. Kubicki, H. Yang, D. Sawada, H. O’Neill, D. Oehme and D. Cosgrove, Sci. Rep., 2018, 8, 13983.

---

Footnote

† Electronic supplementary information (ESI) available: Particle dimensions, demonstration of the reproducibility of titration experiments, comparison of titration curves of dispersed CNC and modified pulp, illustration of linearization analyses of dissolved and immobilized phosphates, estimation of the dissociation constant of disodium hydrogen phosphate, quantification of carbonate contaminants, ionic composition of the simulated experiments, Appendix A, Appendix B, Appendix C, Appendix D, symbols index, MATLAB scripts for the simulated conductometric titration curves. See DOI: https://doi.org/10.1039/d5cp01616a

---