Back Accumulation of Diffusive Gradients in Thin-Films Devices with a Stack of Resin Discs To Assess Availability of Metal Cations to Biota in Natural Waters

Determining species, concentrations, and physicochemical parameters in natural waters is key to improve our understanding of the functioning of these ecosystems. Diffusive Gradients in Thin-films (DGT) devices with different thicknesses of the resin or of the diffusive disc can be used to collect independent information on relevant parameters. In particular, DGT devices with a stack of two resin discs offer a simple way to determine dissociation rate constants of metal complexes from the accumulation of the target metal in the back resin disc. In this work, simple approximate expressions for the determination of the dissociation rate constant are reported and applied to a model Ni nitrilotriacetic complex as well as to Zn complexes in the Mediterranean Osor stream. Once the physicochemical parameters are known, one can plot the labile fraction of the metal complexes in terms of the thickness of the diffusion domain. These plots reveal a strong dependence on the nature of complexes as well as on the characteristics of the diffusion domain, and they are of high interest as predictors of availability to biota whose uptake is limited by diffusion.


Homogeneous and heterogeneous resin models
During the preparation of the Chelex resin discs, while the resin gels are becoming solid, the settling of the beads generates a heterogeneous distribution of the beads in the resin disc which concentrates the binding sites close to the bottom surface. We have used an optic microscope to show the heterogeneous distribution of the resin beads from the profile. Figure S1 Photograph with an optic microscope of a Chelex resin. At the bottom of the picture we can see round shapes that correspond to Chelex 100 beads. On the bottom left of the picture there are different beads at different planes and that is why they are seen more blurred.
Therefore, it is convenient to develop a model that takes into account this heterogeneous distribution of the resin beads in the binding disc of a DGT device.
Let us consider a DGT device with a stack of two resin discs, whose resin beads have settled just in the front part of each resin disc (See Figure S2A). In order to assess for the relevance of the settling, we will use the homogeneous model as reference. In contrast to the heterogeneous model, in the homogeneous model, the binding sites are uniformly distributed within the whole resin disc (See Figure S2B). S4 Figure S2 is a schematic representation of the stack of the two resin discs. The x-axis is perpendicular to the interface between the resin disc and the diffuse disc. The resins are in the range 0<x< r . S5 Figure S2 Domains considered in a heterogeneous distribution (panel A) and a homogeneous distribution of the resin beads (panel B) in a stack of two resin discs (front resin,  r /2<x<  r and back resin <x<  r /2). Orange background indicates the presence of resin sites, while they are absent in the blue background shading.  g is the aggregate thickness of the diffusive disc plus that of the Diffusive Boundary Layer (DBL).  r is the aggregate thickness of the two resin discs.
In the next section, an analytical solution for the diffusion-reaction equations involving heterogeneous resin discs is reported.

Derivation of a solution for cM and cML for a heterogeneous resin of a DGT device
We have developed a model that solves the diffusion-reaction equations for the boundary conditions of a DGT device considering a heterogeneous distribution of the beads of the resin in the resin discs.
Let us consider the simple case of a free metal (M) that reacts with a ligand (L) to produce a complex (ML). The equilibrium constant of this reaction process is defined as: We assume excess of ligand conditions, which imply that the concentration of L (cL) is much larger than the concentration of M (cM), so that it stays approximately constant everywhere in the DGT spatial domain. Therefore, a conditional equilibrium constant can be defined as: Along the heterogeneous distribution presented in Figure S2A, the following assumptions are going to be considered to solve the diffusion-reaction equations: • The resin domain is r 0 x  , but only in the subdomains • The association of M with the resin beads is considered fast and strong, so perfect sink conditions for M are applied. •

Solution
In the gel layer ( r g r x    +   ), which, according to the Figure S2A is referred as region 1, the transport equations are: An equivalent system of equations that is uncoupled can be found by searching linear combinations of eqs (S13) and (S14). By addition, we obtain: After solving a homogeneous second order differential equation, the general solution for eq (S15) is: . Using eqs (S13) and (S14), after some rearrangements we obtain: The general solution for eq (S17) is . Q and P can be redefined in terms of a new couple of constants A' and B' in order to obtain simplified expressions. In the case of the first region, these parameters are set equal to: Then ϕ in the first region has the following expression: Now, to obtain a particular solution for the free metal concentration and complex concentration, we need to apply the boundary conditions (from the first region) to find the particular solutions for eq (S16) and (S21). These conditions correspond to: where the superscript "r" indicates value of a concentration at x= r . These values ( ) r r M ML and cc , will be found later on (e.g. see eq (S63)) Once applied the boundary conditions (S22)-(S25), we obtain the solution for the integration constants: With the integration constants, we can obtain the concentration profiles in region 1 for cM and cML as: , is the ratio of the diffusion coefficient of the metal complex and the free metal.
It is interesting to calculate r ML x c x  =   for the boundary condition involving the flux. From eq (S29), we obtain: This expression is of the interest to apply a boundary condition involving the flux at the right side of the resin-gel interface (eq (S30)) which will allow one to find an expression for r ML c in terms of the physical parameters.
In region 2 ( r r 3 4 x    ), cM=0. The system of differential equations (S13)-(S14) is then reduced to: A general solution for eq (S31) is : Now, to obtain a particular solution for eq (S31), we need to apply the boundary conditions (S25) together with where the superscript "3r/4" indicates that the value of the concentration is taken at x=3 r /4.

S11
The solutions for the integration constants are: Let us, now, consider the region 3 ( rr 3 42 x   ) where eq (S15) and eq (S14) have to be solved. The general solutions for cM and cML are eqs (S16) and (S21).
In region 3, we have to apply eq (S33) and the following boundary conditions for the free metal and complex: With the integration constants, we can write the concentration profiles in region 3 for cM and cML in terms of the boundary values (e.g. r 2 ML c ) and the physicochemical parameters as: From these expressions, later on, we are going to find 3 r 4 ML c and r 2 ML c in terms of the physicochemical parameters.
In region 4, ( rr 24 x   ), the same differential equation as in region 2 has to be solved, (eq (S31)) whose general solution is: In this region 4, the boundary conditions to obtain the particular solution of the differential equation are (S39)   sinh 4 From these expressions we will find an expression for r 2 ML c and r 4 ML c in terms of the physicochemical parameters.
In region 5, (0 <x< δ r /4), the same system of uncoupled differential equations as in region 1 and 3 has to be solved (eqs (S15) and (S17). Their general solutions are eq (S16) and an equation similar to eq (S21) : Now to obtain a particular solution for the free metal concentration and complex concentration in region 5, we need to apply the boundary conditions.
The solutions for the integration constants are: With the integration constants, we obtain the concentration profiles for cML and cM in region 5 as: ( In order to obtain 3r 4 ML c , we apply: from which, using eqs (S37)  c in a compact way, we define the variables a, b and z :

Deriving the back percentage B for a heterogeneous resin
To derive an expression to determine the percentage of back accumulation, B , one needs to calculate the ratio of accumulated moles of metal from x=0 to x= r /2 over the accumulated moles of metal from x=0 to x=  r .
In order to derive this expression, first, one needs to add the differential equations for the free metal and complex (eq (S13) and eq (S14)). Second, one needs to integrate twice the resulting differential equation (S15). After obtaining the general solution (S16) and applying the boundary conditions (S22)-(S25), we obtain: A similar procedure leads to obtain the flux at x= δ r/2 , but one needs to apply the boundary conditions of the region 3 ( eqs (S33), (S38) and (S39)).

Deriving B for a homogeneous resin
As above, to derive an expression for B in a homogenous resin, one needs to calculate the ratio of accumulated moles of metal from x=0 to x=δ r /2 over the accumulated moles of metal from x=0 to x= δ r .
In steady-state and ligand excess conditions, the accumulation when a metal (M) reacts with a ligand (L) to produce a complex (ML) can be written as: With the same procedure used in the heterogenous resin, but using cM =0 for the whole resin domain 0 <x< δ r , the flux at x=δ r can be written as, The main difference in the flux calculation respect to the heterogeneous resin is found in the calculation of the flux of complex at x=δ r /2 . Eq (S31) describes the diffusion-reaction equation for ML in a region with resin beads. This eq holds now for the complex in region 2 of the homogenous resin. The general solution of Eq (S31) can be written as: A particular solution can be found by applying the following boundary conditions:

(S85)
The result for the integration constants A2 and B2 after applying the boundary conditions is :

Determination of kd from B
Sections 5.1 and 5.2 outline the determination of B from the limiting expression (8) of the main manuscript, from eq (S89) or from eq (S78). A short discussion on the validity conditions for using each one of these expression follows. Complementary to these Sections, an additional Excel file is included in the SI ready to estimate kd value using to all these equations once input data from the user is introduced. The most accurate value of kd appears in a green cell (See Section 9 in this file).

Algorithm to obtain kd using the complete expression for B corresponding to a homogeneous resins
The use of the limiting expression eq (8) of the main manuscript, when applicable, allows a straightforward determination of kd, while, when eq (S89) (eq (6) of the main manuscript) needs to be used, an iterative procedure can be applied.
In order to have a simple criterion for the plausibility of the application of eq (8) On the other hand, due to steady state, one can assume that the flux of arriving ML at x= r is also the amount of ML that is dissociating along an effective thickness of the complex penetration into the resin disc given by The fulfillment of eq (S93) allows one to use the limiting expression (8) for the determination of kd. Since both, m and λML depend on kd, we cannot verify eq (S93) until a provisional kdvalue was obtained from eq (8). The fulfilment of eq (S93) with this provisional value ensures that the kd obtained is a good approximation for the dissociation rate constant of the complex.
Conversely, if eq (S93) is not satisfied, a better approximation for kd should be sought using eq (S89) instead of (8).

S24
A simple, but effective, algorithm to calculate kd using eq (S89) is explained as a flowchart in Figure S3. The algorithm follows an iterative procedure that allows to obtain a value for kd from an experimental B (in our work BExp=0.29 in the case of stream Osor). The algorithm searches the kd result that satisfies the Bexp for a given tolerance (Test=1×10 -4 ) within a span of kd values ( 1×10 -1 s -1 to 1×10 -6 s -1 ). This iterative procedure needs an educated guess as an initial value for kd which should be taken as a kd=1×10 -6 s -1 or kd=1×10 -1 s -1 . This educated guess can be assessed depending on the lability degree (ξ) of the complex. For ξ close to 0.7, kd should be close to 1×10 -2 s -1 (as it is expected for a quite labile complex), but for cases where ξ is close to 0.3, kd will be close to 1×10 -4 s -1 (as the behaviour of the complex should be quite inert). Therefore, in the first case a good educated guess would be 1×10 -1 s -1 while for the second case a good educated guess would be 1×10 -6 s -1 . Moreover, for the first case in this algorithm it has to be chosen a negative increment (e.g increment=-1×10 -6 s -1 ), so the kd value decreases and we search for a solution for partially labile complexes. While in the second case a positive increment is needed so the kd increases and then we start searching for a solution of kd for rather inert complexes. Solutions for kd either for very labile or very inert complexes are not feasible with these method as transport is limited by diffusion and the kinetic reaction process does not have an impact on the accumulation.
Alternatively, a value for kd can be obtained from high-level coding programs. For example, Mathematica has a subroutine called Findroot that helps to find a solution for a given equation or system of equations.
With the solution, we need to check whether the heterogeneity of the resin beads distribution is relevant for the determination of kd. To do so, we need to check that Bexp<Bmax, where Bmax is the maximum back percentage that can be reached in a stack of two homogeneous resin discs.
When Bexp is close to Bmax, an accurate computation of kd requires the use of the heterogeneous resin model, by means of eq (S78), as outlined in the following Section.

Algorithm to obtain kd for a heterogeneous resin
According to fig (2) in the main manuscript, eq (S89) for the homogeneous model plotted against kd, has a maximum, that can be labelled as Bmax . Bmax is then, the maximum back percentage that can be reached by a stack of two homogeneous resin discs. Section "Derivation of maximum B for a homogeneous model" gives some details on the calculation of Bmax.
When Bexp is close or greater than Bmax, eq (S78) is needed to recover accurate kd values. In this task, we can use essentially the same algorithm as in section "Derivation of a solution for cM and cML for a heterogeneous resin of a DGT device ". S26 Figure S4 Algorithm applied to determine kd from a B measured in Osor stream for a heterogeneous resin.
In the flowchart shown in Figure S4 Once they are written down in the program, B can be computed as prescribed in eq (S78).
Once again, notice that this procedure is the algorithm that has been used in this work, nevertheless other programs can be applied to iteratively solve the equation and find a solution for kd that satisfies the Bexp.  Figure S5 depicts the concentration profiles of free metal and complex in a stack of two heterogeneous resin discs (as illustrated in Figure S2 ). This figure is quite useful to understand that the percentage of back accumulation, when using heterogeneous resin discs, is higher than the back accumulation when the resin discs are homogeneous. Dissociation of the complex in the whole volume 0 < x < 0.4 mm leads to back accumulation when the resins are homogeneous.
Instead, dissociation in the volume 0 < x <0.5 mm leads to back accumulation when the resins S29 are heterogeneous. Notice in Figure S5 that in the volume 0<x<0.2 mm where there is no resin beads, a profile of free metal resulting from the complex dissociation appears. The decreasing cM concentration as x increases indicates that the free metal produced in this volume diffuses towards increasing x values and gets bound to the back resin at x=0.2 mm. Likewise, a free metal concentration appears in the volume 0.4 mm < x <0.6 mm. The profile of the free metal (continuous blue line in Figure S5) indicates that the metal produced between 0.4 mm < x< 0.5 mm diffuses also towards the back resin disc and binds to it (at x=0.4 mm) increasing the back accumulation with respect to what takes place when the resin discs are homogeneous.

Derivation of maximum B for a homogeneous model
Back accumulation in a DGT device is negligible for a fully labile system as well as for a completely inert one. In both cases there is no dissociation of complexes in the resin domain, which is the phenomenon responsible for the back accumulation assuming that the resin beads are homogeneously distributed in the resin disc with a concentration enough to bind all the metals that diffuse in this domain (i.e. no saturation effects).
The accumulation in the back resin disc depends on the physicochemical parameters of the metal species (DM, DML, ', kd ), as well on the characteristics of a DGT device (δ g and δ r ). In this Section, we are going to derive a condition to determine the dissociation rate constant of The condition of maximum requires that this derivative is equal 0. Unfortunately, there is not an explicit solution of equation (S94) for kd. Therefore, the Bolzano method has been applied to find the kd value (kd,max) which produces a 0 in eq (S94), where B(kd,max)= Bmax .The initial points that define the starting section are selected as kd,a=0.1 s -1 and kd,b=1×10 -9 s -1 . To calculate a middle point kd,c between kd,a and kd,b , we have used log(kd,c)=(log(kd,a)+log(kd,b))/2. The solution for kd is accepted once the condition (kd,akd,b)/2<tolerance is satisfied, where tolerance=1×10 -9 s -1 .
The Bolzano method can be applied iteratively for a range of K' values, so that we can plot a figure of Bmax vs K' We have computed this relationship for the following parameters δ g =1.079×10 -3 m, δ r =8×10 -4 m, DML=DM=4.20×10 -10 m 2 s -1 and the result is shown in Figure   S6 . Figure S6 : Impact of ε ' on Bmax, the dashed black line shows the behaviour of Bmax (computed with the homogeneous model) against log (ε ') for the parameters of δ g =1.079x10 -3 m, δ r =8x10 -4 m, DM=4.20x10 -10 m s -2 and ε=0.1. The solid green line is calculated with the same parameters as the dashed black line but using ε=1.

S31
This figure is of high interest to assess whether the homogeneous model of the resin disc is accurate for determining kd or the heterogeneous model of the resin disc is required. As seen in Fig. 1 of the main text, when the percentage of back accumulation (Bexp) approaches the maximum value of the homogeneous model (Bmax), an accurate value of kd will require the use of eq (S78). Figure S6 allows one to check whether Bexp is close to the maximum B (Bmax ) for a given K' and standard δ g =1.079×10 -3 m and δ r =8×10 -4 m values. Figure S6 also shows that a change of ε from 1 (solid green line) to 0.1 (dashed black line) leads to just small S32 discrepancies in Bmax , when ε ' is smaller than 1. Thus, this plot can be used for the typical values in standard DGT devices.
Just as an example, we can consider the case of Zn in the Osor stream. We have K' 4 70 and =1, so log( K' 67. From Figure S6 , we know that the maximum value that can be achieved for B is near to 0.22. As in our case we have Bexp=0.29 , it is not possible to find a solution for kd using the homogeneous model. Then, instead of applying eq (S89) we will need to apply eq (S78). In the case that Bexp<Bmax , then we can find a solution for kd applying the homogeneous model (eq(S89)) and from our computations we know -for these data-that the kd retrieved by the heterogeneous model will be just 0.3 log units greater than the result retrieved by the homogeneous model which corresponds to a kd value 2 times greater.  If this message does not appear in your screen, please look at internet for how to enable macros for Excel taking into account your Excel version and your operative system.

Structure of the Excel file
The Excel file contains 3 sheets: "Input" "Intermediary calculations" and " utput"

S35
In the Input sheet, you can find the list of nput values that are needed to calculate "kd" By default, the list of entries is filled with the values from the measurements in natural waters in the Osor stream: Figure   he sheet " utput" displa s the odel and the result for kd that is obtained. The correct solution is highlighted in green.