Ferric Ion Diffusion for MOF-Polymer Composite with Internal Boundary Sinks

Simple and economical ferric ion detection is necessary in many industries. An europium-based metal organic framework has selective sensing properties for solutions containing ferric ions and shows promise as a key component in a new sensor. We study an idealised sensor that consists of metal organic framework (MOF) crystals placed on a polymer surface. A two-dimensional diffusion model is used to predict the movement of ferric ions through the solution and polymer, and the ferric ion association to a MOF crystal at the boundary between the different media. A simplified one-dimensional model identifies the choice of appropriate values for the dimensionless parameters required to optimise the time for a MOF crystal to reach steady state. The model predicts that a large non-dimensional diffusion coefficient and an effective association with a small effective flux will reduce the time to steady-state. The effective dissociation is the most significant parameter to aid the estimation of the ferric ion concentration. This paper provides some theoretical insight for material scientists to optimise the design of a new ferric ion sensor.


Introduction
There are many situations in which it is important to monitor the concentration of ferric ions (trivalent iron cation or iron(III)) in a solution. For example, in environmental contexts, a high concentration of ferric ions promotes bacterial and algal growth, which can lead to the death of aquatic animals and plants [1,2]. In the mining industry, the concentration of ferric ions can affect copper and gold yield during mineral leaching, with both high and low concentrations producing adverse effects [3,4]. Other areas in which a knowledge of the ferric ion concentration is important include health and drinking water quality assurance [5,6].
This research is part of a larger project on lean mineral processing wherein a new sensor is being developed to detect ferric ions in a time-, cost-, and energy-efficient manner. Current ferric ion sensing methods are complex, use expensive equipment, suffer from interference from other ions, and are not adaptable to real time monitoring [7]. In addition, ferric ion concentrations are often determined indirectly [7]. As a consequence, a new ferric ion sensor that is stable, cost-effective, energy-efficient, and easy to use is desirable and would impact many sectors. In this paper, we present an idealised model of a ferric ion sensor. We propose a diffusion model of ferric ions through two different media with a sink located at the boundary separating the media.
Xu et al. [8] recently established the effectiveness of a europium-based metal organic framework (EuMOF) for ferric ion sensing. The crystalline EuMOF was suspended in an aqueous solution and tested against various concentrations of metal ions for changes in luminescence. Most MOF ferric ion sensors are intensity-based "turn-off" sensors where the intensity of light emitted by the sensor diminishes in the presence of ferric ions. However, the EuMOF [8] is bimodal, so changes in the emission ratio of two frequency peaks are linearly proportional to the ferric ion concentration, while the intensity of the emission is not as important. Metal ions such as Fe 2+ , Ag 2+ , Ca 2+ , and Zn 2+ have little effect on the emission ratio, and consequently, the ferric ion selectivity of the EuMOF has promising sensing capabilities.
In [9] we investigated the van der Waals interaction energies and forces between a hydrated ferric ion and an EuMOF crystal pore. The Coulombic forces were not considered since ferric ions exist in highly acidic solutions, and the electrostatic forces are negligible [9]. The findings suggest that the hydrated ferric ion is attracted to the pore but does not enter due to steric interactions. This is advantageous from a practical point of view because the hydrated ferric ion can be "washed" away and the sensor can be reused.
In this paper, we investigate how ferric ions diffuse through a solution and a polymer layer and interact with a MOF crystal located on the boundary between the solution and the polymer. One-dimensional diffusion through multiple layers has been studied before. Hickson et al. used a numerical approach to solving diffusion through multiple layers while considering various matching conditions at the boundary [10]. Carr and Turner [11] developed a semi-analytical method to address the complexities that arise for a large number of layers. An analytical solution exists for the one-dimensional problem with two layers (or m layers) for slabs, cylinders, and spheres [12].
A proposed new ferric ion sensor consists of a thin film of polymer coating a cut-away section of optic fibre and embedded with EuMOF crystals. In this scenario, the sensor will have MOF crystals distributed throughout the polymer composite and over its surface. In our model, we place an EuMOF crystal at the boundary between a solution and a polymer, and we investigate the changing concentration of ferric ions bound to the crystal. The paper provides important insight into ferric ion behaviour for a proposed sensor that is still under development and for which there is no readily available experimental data. Squires et al. [13] explore how the physical attributes of a sensor affect the flux onto the sensor, the rate at which substance of interest binds to the sensor, and the approximate times scales for the system to reach equilibrium. They present simple rules to equip the reader with a basic understanding of the environment and enable the design of a sensor better suited to the environment. Yariv [14] developed an advection-diffusion-reaction model of analytes binding to a sensor located on a solid surface adjacent to a shear flow of solution, using an equation first proposed in [13] to describe the concentration on the sensor. Here, we use a similar model to represent the association and disassociation of ferric ions to the EuMOF crystals.
In Section 2, we describe the geometry and the physical environment of the ferric ion sensor. We present a mathematical model that describes both the diffusion of ferric ions through an analyte solution and a polymer region, as well as the association of ferric ions to a MOF crystal located at the boundary between the two regions. Section 3 discusses outcomes of the model and the role and importance of each of the parameters.

Mathematical Modelling and Assumptions
One possible design for a ferric ion sensor is to suspend EuMOF crystals in a polymer composite to create a thin film on a cut away section of an optic fibre, as described by [15]. The ferric ions bind to the MOF pores on the crystal surface. A light pulse is sent through the MOF-polymer-coated optic fibre and a detector measures changes in luminosity. Changes in luminosity are indicative of the concentration of ferric ions in the solution. Here, we are particularly interested in the diffusion of the ferric ions through the mixed media and their association with the MOF crystals in order to gain some insight into which physical parameters might be most important. As such, we study an idealised version of the device depicted in Figure 2. We model only the process of diffusion in the analyte solution, the association of ferric ions to a MOF crystal located at the interface between the polymer composite and the solution, and diffusion in the polymer composite. We position one MOF crystal at the interface between the solution and the polymer matrix. We do not consider the association of ferric ions with MOF crystals inside the polymer composite.
The MOF-polymer composite is exposed to a solution containing ferric ions, where the movement of ferric ions is governed by Brownian motion. We assume that the diffusivity of the ferric ions is constant through the solution, denoted by D 1 . Since, in practice, the polymer is kept hydrated and the ferric ions do not enter the MOF pores, we assume that the diffusivity in the polymer is also constant, D 2 < D 1 , and the Vrentas-Duda theory is not applicable. As a consequence, the movement of ferric ions in both the solution and the polymer is assumed to be governed by the conventional diffusion equation with distinct diffusivities, where C(T, X) is the concentration of ferric ions, and the origin is located at the centre of the MOF strip so that the X 1 -axis coincides with the solution-polymer boundary. Here, ∇ 2 is the usual Laplacian in two dimensions. The ferric ion concentration profile in the analyte solution is given by C(T, X) for X 2 ≥ 0, and the concentration in the polymer is given by C(T, X) for X 2 < 0. The time since the sensor's exposure to the solution is given by T. . Two-dimensional system: X 1 -axis is solution-polymer boundary, X 2 > 0 is distance into solution, and X 2 < 0 is into polymer composite.

Ferric ion diffusion through solution
To ensure that continuity is maintained for boundaries shared by two media, we include two additional conditions. The ferric ion concentration at the boundary is assumed to be the same when approaching from above or below the boundary, and the flux out of the solution is equivalent to the flux into the polymer [10], Far from the sensor in the direction of the analyte solution, we assume a far field boundary condition, where the ferric ion concentration remains constant and equal to that of the bulk solution C 0 . For numerical purposes, we set the location of this condition to be at X 2 = L, so that where L is assumed to be a large distance away from the solution-polymer boundary. We also assume no flux at a horizontal distance L away from the MOF centre (symmetry) and at the bottom boundary of the polymer composite, The flux condition at the solution-MOF interface is given by where B(T, X 1 ) is the occupancy of the MOF pores by the ferric ions, k on is the association rate compared to the ferric ion concentration, k off is the dissociation rate, B 0 is the number of MOF pores available for ferric ion occupation, and a is a constant. Due to steric effects, the ferric ion only associates at the MOF pore's entrance and does not diffuse into the MOF pore. The first term captures the association of ferric ions to the MOF crystal, and this term vanishes when the MOF crystal is completely occupied. The second term captures the dissociation of ferric ions from the MOF pore. We assume that ferric ions only associate and dissociate from the MOF pores in the X 2 -direction and impose a no-flux boundary condition to the sides of the MOF crystal, The occupancy of the MOF pores is governed by the pseudo-ordinary differential equation (ODE), In addition, since B 0 is the number of MOF pores available, B(T, X 1 ) ∈ [0, B 0 ] and the total ferric ion occupancy in the MOF crystal is, The assumed initial conditions are

Dimensionless Two-Dimensional Model
We non-dimensionalise distance with the half length of the MOF strip, , the ferric ion concentration with the bulk concentration, C 0 , and the MOF binding site occupancy with the maximum occupancy, B 0 . There are two options for the time scale: a diffusive time-scale, where τ = 2 /D 1 , or an association time-scale, where τ = 1/C 0 k on . At this stage, we non-dimensionalise time with τ without specifying which time-scale so that the system of equations can be written where D = D 2 /D 1 , s 1 = D 1 τ/ 2 , s 2 = C 0 k on τ and s 3 = k off /C 0 k on are non-dimensional parameters. If the problem is scaled with the diffusive time-scale, then s 1 = 1. If the problem is scaled with the association time-scale, then s 2 = 1.
The continuity conditions at the solution-polymer boundary become Far away from the solution-polymer boundary, the boundary condition becomes while the no-flux boundary conditions at the sides and at the polymer-optic fibre boundary become At the solution-MOF boundary, where s 4 = a B 0 / C 0 is a non-dimensional constant, and the no-flux boundary conditions at the sides of the MOF crystal become The initial conditions and the total MOF pore occupancy are

One-Dimensional Model
To obtain a preliminary indicator of the behaviour of the system and the importance of the various parameters, we further simplify this model. Figure 3 shows a simpler onedimensional system, where the MOF crystal is represented by a red square, and diffusion in the polymer composite region is completely ignored.

Column of solution
Polymer surface The concentration of ferric ions is assumed to be only a function of two variables, C(T, X), and the MOF pore occupancy is assumed to depend only on time, B(T), The corresponding flux condition at the boundary X = 0 is given by and the far field condition is C(T, L) = C 0 , with initial conditions C(0, X) = C 0 and B(0) = 0.

Dimensionless One-Dimensional Model
We non-dimensionalise concentration, occupancy, length, and time in the same way as before, where is a characteristic length representative of the MOF crystal size so that the system of equations becomes where the s i are as previously defined. The non-dimensional flux condition at the boundary x = 0 is while the non-dimensional far field condition at x = L/ becomes and the initial conditions are given by

Results and Discussions
In this section, we discuss the results for the diffusion of ferric ions through the two regions and the ferric ion occupancy in the MOF pores. We vary the dimensionless constants, s i and D, to analyse how the parameters affect ferric ion diffusion and MOF occupancy.
The models are solved numerically, as no analytical solution exists for this system, using purpose written code in MATLAB. First, we examine the dimensionless one-dimensional model where the equation for b(t) is solved using a Runga-Kutta method, and the diffusion equation for c(t, x) is solved with an implicit Crank-Nicolson scheme. Second, we examine the dimensionless two-dimensional model where the equation for b(t, x 1 ) is solved using a Runga-Kutta method, and the diffusion equation for c(t, x) is solved with a forward-time-centred-space finite difference scheme.
Further details about the numerical schemes for the dimensionless one-dimensional and two-dimensional models are presented in Appendices A and B. In addition, the authors have tested some extremely simple cases; for example, if the MOF crystal is not present. However, these results are trivial and are not included here.  Figure 4b shows the MOF occupancy. Ferric ions in the solution initially associate rapidly to the MOF crystal, reducing the concentration in the solution near the MOF crystal. Over time, the ferric ion concentration returns to the steady-state far-field condition, as shown in Figure 4a. We note that the MOF crystal's occupancy initially increases rapidly and then slowly increases until steady state is reached, as shown in Figure 4b.

Dimensionless One-Dimensional Model
In the following subsections, we describe the impact of varying parameters. In all cases for the one-dimensional scheme, ∆x = 0.1 and ∆t = 0.005, and we set L/ = 10.   Figure 5 shows the effect of varying the non-dimensional diffusion coefficient (s 1 = D 1 τ/ 2 ) on the ferric ion concentration in the solution at the MOF crystal and the MOF occupancy. A smaller non-dimensional diffusion coefficient causes the concentration of ferric ions in solution to be more significantly reduced near the MOF boundary, and, following this initial reduction, it takes longer for the concentration in the solution to recover towards the steady-state concentration. A larger non-dimensional diffusion coefficient results in less reduction near the boundary because ferric ions can more readily diffuse from the column of solution. MOF occupancy reaches close to steady-state sooner with a high non-dimensional diffusion coefficient.

Varying the Effective Association Parameter, s 2
Increasing the effective association parameter (s 2 = C 0 k on τ) means that ferric ions associate to the MOF crystal faster than they diffuse through the solution. This results in the MOF occupancy reaching the steady state sooner, as shown in Figure 6b. When s 2 is large, the concentration near the solution-polymer boundary is reduced very quickly, followed by a slow increase as ferric ions diffuse from the rest of the region, as shown in Figure 6a.
Upon exposure of the MOF to the solution, the concentration at the boundary varies significantly for different values of s 2 . However, at large times, the ferric ion concentration and the MOF occupancy is largely independent of s 2 .

Varying the Effective Dissociation Parameter, s 3
The effective dissociation parameter (s 3 = k off /C 0 k on ) captures the ferric ions that dissociate from the MOF crystal. We have investigated various values of s 3 , including s 3 = 0, which represents the situation when ferric ions can only associate-see Figure 7.
Increasing the effective dissociation parameter means that the ferric ions dissociate more readily and the non-dimensional steady-state MOF pore occupancy can be calculated from Equation (11), A similar expression for the number of effective bound receptors at equilibrium is given in Squires et al. [13], with a different combination of parameters. (The present authors believe that the equilibrium constant K D in Squires et al. [13] should be defined as K D = k off /k on rather than K D = k on /k off as stated. This also aligns with working in their Box 1). In the steady state, the occupancy of the pores is lower for larger values of the effective dissociation parameter, and the steady state is reached more quickly. In the case when the ferric ions cannot dissociate, that is s 3 = 0, the steady state MOF occupancy is equal to its maximum value (b(t → ∞) = 1), and the concentration of ferric ions at the MOF-solution boundary approaches the bulk concentration at large times.

Varying the Effective Flux Parameter, s 4
An important consequence of the flux condition is the time taken for the MOF crystal to approach steady state. Increasing the effective flux parameter (s 4 = aB 0 / C 0 ) increases the time taken to reach the steady state (see Figure 8b). The increase affects ferric ion concentration at the MOF crystal boundary, as ferric ions quickly associate to the MOF crystal and the ferric ions do not diffuse fast enough for the MOF crystal to reach steady state in a timely manner, as shown in Figure 8a.
The delay to reaching steady state is attributed to ferric ions associating to the MOF crystal faster and quickly depleting the ferric ions in the solution. A comparable mechanism is observed when increasing the effective association parameter, s 2 . The effective association parameter is a MOF characteristic, and an increase results in a shorter time for the MOF crystal to reach steady state (see Figure 6b). The ferric ion concentration at the MOF boundary is largely unchanged after t = 4, as shown in Figure 6a. However, the effective flux parameter is a ferric ion characteristic, and an increase in this parameter results in a longer time for the MOF crystal to reach steady state, as shown in Figure 8b. The ferric ion concentration at the MOF boundary is very sensitive to the increase in the effective flux parameter, as seen by the long depletion time for a larger parameter value in Figure 8a.

Dimensionless Two-Dimensional Model
In this section, we return to the dimensionless two-dimensional model and to the primary aim of the paper, which is to investigate ferric ion diffusion through a solution into a polymer composite with a MOF sink located at the boundary between the two media. Figure 9 shows the ferric ion concentration for the domain −5 ≤ x 1 ≤ 5, −5 ≤ x 2 ≤ 5, where x 2 < 0 is the concentration in the polymer and x 2 ≥ 0 is the concentration in the solution. The MOF crystal is located at |x 1 | ≤ 1 and x 2 = 0. Figure 9a,b show the ferric ion profiles at time t = 0.5 and 5, respectively. In all cases for the numerical twodimensional scheme, ∆x 1 = ∆x 2 = 0.1 and ∆t = 0.0005, and D = 0.5 (unless otherwise stated) and L/ = h/ = 10. Figure 10a shows the ferric ion concentration along the line x 1 = 0, which shows the changes through the mixed media over time. In comparison to the one-dimensional model (Figure 4a), the two-dimensional model requires additional time to reach steady state. This is due to the size of the MOF crystal and the time taken to reach steady state in the two-media system. The MOF occupancy at the centre of the MOF (b(t, 0)) and the total MOF occupancy over the whole MOF crystal (Equation (10)) is shown in Figure 10b. By comparing Figure 10b with Figure 4b, we see that the MOF occupancy as predicted by the two-dimensional model is very similar to that predicted by the one-dimensional model. However, the two-dimensional model includes the effect of the polymer on ferric ion diffusion and greater MOF occupancy potential.

Varying the Non-Dimensional Parameters
Variation of the dimensionless constants has the same effect on the two-dimensional MOF occupancy b(t, 0) profile when compared to the one-dimensional MOF occupancy b(t). The two-dimensional ferric ion concentration at the centre of the MOF crystal behaves in a similar way to that of the one-dimensional ferric ion profile shown in Figure 4a. However, more time is needed for the ferric ion concentration to reach steady state due to a larger MOF crystal surface area. Figure 11 shows the effect of varying the non-dimensional diffusion coefficient (s 1 = 0.1, 1, 5). For all values of s 1 , the concentration in the solution near the MOF surface initially decreases before increasing towards the steady state value. For larger values of s 1 , the decrease is more significant, which is in contrast to the one-dimensional results, as shown in Figure 5.
The changes in the other non-dimensional parameters (s 2 , s 3 , s 4 ) produce similar results to those discussed for the one-dimensional model. Figures showing the details can be found in Appendix C.

Varying the Relative Diffusion Coefficient, D
If ferric ion diffusivity in the solution and the polymer is the same, then the relative diffusion coefficient is unity, D = D 2 /D 1 = 1. In practice, the diffusion coefficient of ferric ions in the solution will exceed that in the polymer and the larger the difference the smaller the relative diffusion coefficient. Figure 12 shows the effects of varying the relative diffusion coefficient.
Decreasing the relative diffusion coefficient means that ferric ions take longer to diffuse in the polymer and away from the solution-polymer boundary. This affords the ferric ions more opportunity to associate to the MOF crystal. Figure 12c shows that the occupancy in the crystal is highest in the case of the smallest relative diffusion coefficient.

Conclusions
In this paper, we have modelled a proposed new ferric ion sensor. We first examined a two-dimensional model to describe the diffusion of the ions through a solution and into a polymer layer, as well as their attachment to an EuMOF crystal located at the boundary of the two media. We then examined a simplified model that includes diffusion through the solution only, together with attachment at the MOF crystal. The solutions to the corresponding non-dimensional models are then analysed to investigate the importance of the various non-dimensional parameters.
The one-dimensional model indicates that at steady state, the occupancy of the MOF crystal is given by Equation (15), and this occupancy depends only on the initial concentration of ferric ions in the solution and the effective dissociation constant, s 3 = k off /C 0 k on . As a consequence, if the steady-state occupancy, B eq = B(t → ∞), can be measured with the new sensor, and the association rate k on , dissociation rate k off , and maximum occupancy B 0 are known, then the concentration of ferric ions in solution can be calculated (in dimensional terms) using In practice, it is most likely to be useful to minimise the time taken to reach steady state so that sensing of multiple samples can be carried out quickly and efficiently. To achieve this, it would be helpful to have a large non-dimensional diffusion coefficient s 1 and effective association constant s 2 with a small effective flux constant s 4 .
While the two-dimensional model is more representative of the proposed sensor and provides more information about the behaviour of the ferric ions in the solution and the polymer, many of the broad behaviours and impacts of changing the system parameters are captured by the one-dimensional model. The relative diffusion coefficient D suggests that the polymer composite should hinder ferric ion diffusion in that region. A final sensor design is likely to include MOF crystals embedded within the polymer matrix rather than on the surface only, so the polymer composite should be designed so that the relative diffusion coefficient is close to unity in order for ferric ions to associate to MOF crystals inside the polymer composite. This situation is not examined here but will be the subject of future work. In this case, a two-dimensional model is likely to provide more insight. Data Availability Statement: Numerical details for this study are presented in the Appendices. Data sharing is not applicable for this article.

Acknowledgments:
The authors wish to gratefully acknowledge Linda Rozenberga for much chemical knowledge and her considerable insight.

Conflicts of Interest:
The authors declare that they do not have any conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

MOF
Metal organic framework EuMOF Europium based metal organic framework ODE Ordinary differential equation Substituting these formulas into the ferric ion diffusion equation gives We collect terms with time step n + 1 on the left and n on the right hand side of the equation so that the general scheme is where α = s 1 ∆t/∆x 2 . The boundary conditions concerning the ferric ion concentration need to be incorporated, and so, the flux at the MOF-polymer boundary (Equation (12)) becomes for i = I, far-field boundary (Equation (13)), c n I = 1, ∀n, The finite difference scheme can be written as a linear system, , A = [4∆xs 2 s 3 s 4 b n+1 /s 1 , 0, · · · , 0, 2α, 1] T , and c n = [c n 0 , c n 1 , c n 2 , · · · , c n I−2 , c n I−1 , c n I ] T , where L and R are real (I + 1) × (I + 1) matrices, and A and c n are (I + 1)-vectors.
The system describing the concentration of ferric ions through the solution can be evaluated as where c(0, x) = c 0 i = 1. Table A1 provides the numerical values for the graphs shown in the main text. Tests were performed for various timesteps and mesh sizes to ensure numerical accuracy.

Appendix B. Algorithm for Dimensionless Two-Dimensional Model
The functions describing the ferric ion concentration c(t, x) and MOF occupancy b(t, x 1 ) are discritised, and the following definitions are used: for i = ir + 1 and j = jm, solution-polymer continuity condition solution-polymer continuity condition (Equation (3)) and no-flux boundary condition (Equation (7)) at x 1 = 1 + and x 2 = 0, c n+1 ir+1 jm = (1 − 2s x − s y (1 + D))c n ir+1 jm + s x (c n ir+2 jm + c n ir+2 jm ) + s y c n ir+1 jm+1 + s y Dc n ir+1 jm−1 , At the solution-MOF-polymer boundary, we use the continuity condition (Equation (3)), to find  . Table A2 provides the numerical values for the graphs shown in the main text. Tests were performed for various timesteps and mesh sizes to ensure numerical accuracy. Unless otherwise stated, D takes the value found in Table A2.

Appendix C. Results of Varying Parameters in the Two-Dimensional Model
This section contains the figures that show the effects of variation of the remaining non-dimensional parameters for the two-dimensional model.