Reaction-Multi Diffusion Model for Nutrient Release and Autocatalytic Degradation of PLA-Coated Controlled-Release Fertilizer

A mathematical model for the reaction-diffusion equation is developed to describe the nutrient release profiles and degradation of poly(lactic acid) (PLA)-coated controlled-release fertilizer. A multi-diffusion model that consists of coupled partial differential equations is used to study the diffusion and chemical reaction (autocatalytic degradation) simultaneously. The model is solved using an analytical-numerical method. Firstly, the model equation is transformed using the Laplace transformation as the Laplace transform cannot be inverted analytically. Numerical inversion of the Laplace transform is used by employing the Zakian method. The solution is useful in predicting the nutrient release profiles at various diffusivity, concentration of extraction medium, and reaction rates. It also helps in explaining the transformation of autocatalytic concentration in the coating material for various reaction rates, times of reaction, and reaction-multi diffusion. The solution is also applicable to the other biodegradable polymer-coated controlled-release fertilizers.


Introduction
Controlled-release fertilizer (CRF) contains many important sources of nutrients, such as nitrogen, phosphorus, and potassium, that are crucial to plants [1,2]. The CRF granules are coated with polymers [3], organic polymers, biopolymers, natural biomolecule materials, and nanocomposites [2]. The coating materials play an important role in nutrient release profiles by providing the obstacles to the mixing of nutrients with water and help in the water retention of the granules [4]. Biodegradable materials have shown good results in maintaining the optimal control release rates and are environmentally-friendly, as well [5]. Some of the biodegradable material [6], such as starch [6], polyvinyl alcohol [7], 3,5-pyridinedicarboxylic acid N-oxide [8], and poly(lactic acid) [9] and nanocomposites and nanomaterials [10,11] are used as a coating of CRF granules. PLA, a natural synthetic polymer available in abundance, has shown good results in biodegradability, controlled-release rates, and is used extensively in drug delivery and implants [12][13][14][15].
Biodegradation of biopolymers is a complicated process as it occurs through scission of the main chains and side chains of macromolecules. Biodegradation is induced by many factors, such as thermal activation, hydrolysis, biological activity or enzymes, oxidation, and radiolysis [16]. PLA is an aliphatic polyester and is prone to chemical hydrolysis of ester backbones in the aqueous environment.
This makes it difficult to design and optimize the nutrient release applications. There have been many studies with regard to PLA modeling, but most of the papers dealt with the drug delivery applications in pharmaceuticals [15]. There is limited research that deals with the modeling of PLA-coated CRF. For the case of CRF, the nutrient release rate from granules with different size and diffusivity depends mainly on release regime. Hence, the study of nutrient release should necessarily be undertaken in conjunction with the knowledge of polymer degradation.
Modeling of CRF in the single diffusion mechanism was considered by Jarrel et al. [21][22][23]. They discovered that the nutrient release profile was coming from the sulfur-coated urea. Later, Al-Zahrani took into account the single solute diffusion mechanism to predict the nutrient release from CRF [24]. Shaviv et al. [25,26] considered non-Fickian diffusion to model the nutrient release from CRF. In these studies, the nutrient release was divided into three stages; that is, lag period, linear release, and extended release stages for singular granule and for the population of the granules. Lu et al. [27] developed an explicit mathematical analysis of release from coated particle using multi-diffusion phenomenon. Trinh et al. [28][29][30][31] used the multi-diffusion model by employing the finite element method to investigate the nutrient release in the decay period (extended release stage). They also studied the effect of imperfect coating on nutrient release in both soil and water.
For the case of modeling of autocatalytic degradation reaction of PLA in the literature, the models used the homogeneous degradation before the starting of the erosion [32,33]. A PLA-based reaction-diffusion model has been proposed by Wang et al. [34] for autocatalytic degradation for polymer implants in the form of plates and cylinders. The model can be used to other geometries by solving numerically using finite element analysis for the degradation map. The other models of polymer degradation [17,35,36] assumed that hydrolytic degradation follows the pseudo-first-order kinetics to include the autocatalytic effect.
To the best of our knowledge with regard to CRF modeling, there is no model available in the literature, which can predict the nutrient release profile for the reaction-multi diffusion system. The purpose of the current study is to investigate the effect of reaction rates coupled with multi-diffusion on the nutrient release profiles and to obtain insight of autocatalytic degradation of the polymer coating. For this purpose, an analytical-numerical solution was derived from a reaction-multi diffusion model by employing the Laplace and Zakian method. This model can also be extended to other studies of CRF granules using biodegradable polymers as a coating material.

Mathematical Model
A spherical, coated fertilizer granule has been considered as portrayed in the schematic diagram ( Figure 1). It has a radius (a), uniform coating thicknesses (p), and a radius of the coated granule (b). The initial nutrient concentration in the granule (C 0 ) which is, in turn, less than nutrient saturation concentration (C s ). The nutrient concentration inside the granule, coating film, and extraction medium are represented by C m (r, t), C f (r, t), and C e , respectively. The diffusion coefficient of the granule and coating thickness are denoted as D m and D f , respectively.
When the PLA-coated granules are subjected to dissolution in wet soil, the water molecules penetrate through the coating membrane, and dissolve the nutrient which gives rise to the osmotic pressure inside the coating shell. The nutrient then diffuses from the dissolved core, through the coating membrane, and into the external environment. This chemical phenomenon is represented by the mathematical equations in Equations (2) and (3). The basic equations of multi-diffusion model along with reaction in a sphere are given as follows [37]: The initial conditions are [27]: (0) = 0 (4c) and boundary conditions are [27]: and are the constant values used in the partition. is the volume of extraction medium in which the nutrient dissolves. In Equation (3), is given as the net generation of species per volume. The chemical reaction for interest is autocatalysis of PLA degradation. For the PLA autocatalysis, when the concentration of the ester bonds and water is constant, it follows the pseudofirst order rate law [38]. Siparsky et al. [17] conducted an experimental study on PLA hydrolysis in aqueous medium and concluded that the PLA hydrolytic degradation shows good agreement with When the PLA-coated granules are subjected to dissolution in wet soil, the water molecules penetrate through the coating membrane, and dissolve the nutrient which gives rise to the osmotic pressure inside the coating shell. The nutrient then diffuses from the dissolved core, through the coating membrane, and into the external environment. This chemical phenomenon is represented by the mathematical equations in Equations (2) and (3).
The basic equations of multi-diffusion model along with reaction in a sphere are given as follows [37]: The initial conditions are [27]: C e (0) = 0 (4c) and boundary conditions are [27]: ).
(5e) K a and K b are the constant values used in the partition. V e is the volume of extraction medium in which the nutrient dissolves. In Equation (3), R v is given as the net generation of species per volume. The chemical reaction for interest is autocatalysis of PLA degradation. For the PLA autocatalysis, when the concentration of the ester bonds and water is constant, it follows the pseudo-first order rate law [38]. Siparsky et al. [17] conducted an experimental study on PLA hydrolysis in aqueous medium and concluded that the PLA hydrolytic degradation shows good agreement with pseudo-first-order Polymers 2017, 9, 111 4 of 13 chemical kinetics. However, once the reaction is taken into consideration, the concentrations vary depending upon time and reaction rate. Hence: For the purpose of generality, the above problem is solved using the non-dimensional technique as defined in Equations (7a)-(7c): 3 3 , The conservation equation is given by Equations (8) and (9): The initial and boundary conditions are also transformed as represented in Equations (10a)-(10h): and

Analytical Solution of the Model
To solve the above set of equations represented in Equations (8) and (9) with the help of initial and boundary conditions given from Equations (10a)-(10h). The proposed solution is given by using the Laplace transform, which is defined as: Polymers 2017, 9, 111 5 of 13 By using the above definition of Laplace transform, Equations (8) and (9) can be transformed as: Boundary conditions are also transformed and written as: The solution of Equations (12) and (13) is given by the following results. The detailed derivation of the solution is not shown here due to lengthy steps. In the current paper, only significant steps are displayed. Thus, Equations (12) and (13) can be transformed into Equations (15) and (16), as shown below: where A 1 , A 2 , B 1 , and B 2 are the constants to be determined using the initial and boundary conditions from Equations (14a)-(14e). By using Equation (14a) in Equation (15), we get B 1 = 0. Then, Equation (15) can be reduced to Equation (17): To find the constants A 1 , A 2 , and B 2 , the detailed procedure is presented from Equations (A1)-(A7) as presented in the Appendix A.

Numerical Inversion of Laplace Transform
The above equations have to be inversed to obtain the solution. However, as these equations are complex and do not have a pole at s = 0, it is very difficult to obtain the inverse Laplace using the residue theorem [39]. The above problems can be solved by using numerical inversion of the Laplace transformation given by Zakian [40]. Zakian has developed an explicit formula for numerical inversion using a weighted function to approximate the time domain function.
where K j and α j are the constant in complex conjugate pairs. N is the number of terms used in the summation. f (t) approaches f (t) only when N reaches ∞. The previous study has shown that the truncated error is negligible when N = 5 [41]. The constant values for K j and α j are shown in Table 1. The dimensionless cumulative release of g(τ) is equal to M t /M ∞ , which is equal to θ e (τ)/θ e (∞). At infinity, the reaction rate of autocatalysis degradation is negligible. Hence, θ e (∞) is calculated as given by Lu et al. [27], which can be represented in Equation (22). The autocatalytic concentration profiles are obtained by solving the Equation (19):

Nutrient Release Profiles
The simulation of reaction-multi diffusion has been carried out for a dimensionless coating thickness l = 1, different volume of extraction medium, V r = 1, 3, 10, K a = K b = 1, and different dimensionless diffusivity, D r = 0.01, 0.1, 1, and 10. The nutrient release profiles have been presented from Figures 2-4. In Figure 2, the release profiles have been portrayed only for multi-diffusion. In this case, the reaction rate constant (k) has been considered as zero. The model equations are presented by Equations (23a) and (23b). These equations can be solved with appropriate initial and boundary conditions from Equations (10a)-(10h). The diffusivity of a material can hold the nutrient from release. To see these effects, various D r are chosen, whereby D r = 0.01, 0.1, 1, and 10, while V r = 10. The dimensionless nutrient release with dimensionless time has been shown in Figure 2. As the diffusivity increases, the nutrient release rate also increases [43]. The time taken for 99% of nutrient release (τ 99 ) at D r = 10 is less than τ 99 when D r = 1. The same profile can be seen when D r = 0.1 and 0.01.  Figure 4. The value of the extraction medium plays an important role. For example, when the value of diffusivity is less than 1 [27], its result can be clearly seen from Figures 3 and 4. It can also be inferred that when the reaction rate k = 0, the 90 is larger for both the cases as explained in Figures 3 and 4 as compared to 90 when k = 0.015. When the reaction rate constant (k) increases, 90 decreases. This indicates more release of nutrients in the shorter duration of time. When the initial stage of nutrient release is taken as 30%, the time taken ( 30 ) is the same for all the possible k values. These phenomena occur because the release of nutrient at the initial stage is taking place when the degradation is negligible. However, it makes a significant effect in the second half ( 40 ) of the nutrient release process as the degradation takes place by creating more pores in the coating, which encourages the fast release of nutrients from the granule. Hence, it has helped to understand the effect of reaction along with diffusion in the nutrient release profiles instead of dependence only on the diffusion coefficient and diffusion phenomenon as considered by Lu et al. [44].   Nutrient release profiles by considering both reaction and multi-diffusion systems have been obtained for different reaction rates, that is, k = 0, 0.015, 0.03, 0.02, 0.04, 0.06, 0.08 day −1 . The simulation was carried out for constant diffusivity, whereby = 1, = 1, and = 10. The results are shown in the Figure 3. For the case of = 0.1 and = 1, the results are shown in Figure 4. The value of the extraction medium plays an important role. For example, when the value of diffusivity is less than 1 [27], its result can be clearly seen from Figures 3 and 4. It can also be inferred that when the reaction rate k = 0, the 90 is larger for both the cases as explained in Figures 3 and 4 as compared to 90 when k = 0.015. When the reaction rate constant (k) increases, 90 decreases. This indicates more release of nutrients in the shorter duration of time. When the initial stage of nutrient release is taken as 30%, the time taken ( 30 ) is the same for all the possible k values. These phenomena occur because the release of nutrient at the initial stage is taking place when the degradation is negligible. However, it makes a significant effect in the second half ( 40 ) of the nutrient release process as the degradation takes place by creating more pores in the coating, which encourages the fast release of nutrients from the granule. Hence, it has helped to understand the effect of reaction along with diffusion in the nutrient release profiles instead of dependence only on the diffusion coefficient and diffusion phenomenon as considered by Lu et al. [44].  Nutrient release profiles by considering both reaction and multi-diffusion systems have been obtained for different reaction rates, that is, k = 0, 0.015, 0.03, 0.02, 0.04, 0.06, 0.08 day −1 . The simulation was carried out for constant diffusivity, whereby D r = 1, l = 1, and V r = 10. The results are shown in the Figure 3. For the case of D r = 0.1 and V r = 1, the results are shown in Figure 4. The value of the extraction medium plays an important role. For example, when the value of diffusivity is less than 1 [27], its result can be clearly seen from Figures 3 and 4. It can also be inferred that when the reaction rate k = 0, the τ 90 is larger for both the cases as explained in Figures 3 and 4 as compared to τ 90 when k = 0.015. When the reaction rate constant (k) increases, τ 90 decreases. This indicates more release of nutrients in the shorter duration of time. When the initial stage of nutrient release is taken as 30%, the time taken (τ 30 ) is the same for all the possible k values. These phenomena occur because the release of nutrient at the initial stage is taking place when the degradation is negligible. However, it makes a significant effect in the second half (τ 40 ) of the nutrient release process as the degradation takes place by creating more pores in the coating, which encourages the fast release of nutrients from the granule. Hence, it has helped to understand the effect of reaction along with diffusion in

Autocatalytic Concentration Profiles
The analytical expression gives the importance of first order reaction rates. This helps to understand whether the PLA degradation and erosion is increasing or decreasing by diffusion. However, the degradation reaction cannot fully complete as the ester bonds in the PLA are converted into the monomers after the hydrolysis reaction. The maximum time required for the hydrolysis reaction to take place can be calculated using Equation (24) [27].
The maximum time required to convert all the ester bond in the PLA to monomers is given below: where 0 is the initial molecular weight and M1 is the average weight of monomers. Piemonte et al. [20] found that the hydrolytic degradation of PLA (k = 0.081 day −1 ) can be calculated at temperatures higher than 140 o C . Since the simulation carried out in the current study considers the PLA degradation in the aqueous environment, which is equal to the room temperature, so, we chose 0.015 day −1 as the initial point of the range of degradation. The autocatalytic concentration for degradation reaction without diffusion is described by [45]: The solution of Equation (25) is given below: where is referred as reaction dominant limit commonly used in autocatalytic degradation of PLA [46]. 0 is the initial concentration of monomers in the PLA coating undergoing hydrolytic degradation. The effect of reaction dominant limit and dimensionless concentration profile is shown in Figure 5.
The dimensionless autocatalytic profiles of the reaction multi-diffusion equation has been obtained with time when = 0.1, = 3 and the space radial distance, = 2 and k = 0.05 day −1 . The initial amount of monomers in the PLA, Ct0 = 1.73 × 10 4 mol/m 3 [19]. As the degradation and diffusion dominate in the initial time until the monomers are converted into oligomers, the dimensionless autocatalytic concentration reduces exponentially. The concentration profiles have been taken at = 2 because after that the concentration gets diffused into the extraction medium. The simulation of dimensionless autocatalytic concentration profiles with dimensionless time are shown in Figure 6.

Autocatalytic Concentration Profiles
The analytical expression gives the importance of first order reaction rates. This helps to understand whether the PLA degradation and erosion is increasing or decreasing by diffusion. However, the degradation reaction cannot fully complete as the ester bonds in the PLA are converted into the monomers after the hydrolysis reaction. The maximum time required for the hydrolysis reaction to take place can be calculated using Equation (24) [27].
The maximum time required to convert all the ester bond in the PLA to monomers is given below: where M n0 is the initial molecular weight and M 1 is the average weight of monomers. Piemonte et al. [20] found that the hydrolytic degradation of PLA (k = 0.081 day −1 ) can be calculated at temperatures higher than 140 • C. Since the simulation carried out in the current study considers the PLA degradation in the aqueous environment, which is equal to the room temperature, so, we chose 0.015 day −1 as the initial point of the range of degradation. The autocatalytic concentration for degradation reaction without diffusion is described by [45]: The solution of Equation (25) is given below: where θ r is referred as reaction dominant limit commonly used in autocatalytic degradation of PLA [46]. θ r0 is the initial concentration of monomers in the PLA coating undergoing hydrolytic degradation. The effect of reaction dominant limit and dimensionless concentration profile is shown in Figure 5.
The dimensionless autocatalytic profiles of the reaction multi-diffusion equation has been obtained with time when D r = 0.1, V r = 3 and the space radial distance, η = 2 and k = 0.05 day −1 . The initial amount of monomers in the PLA, C t0 = 1.73 × 10 4 mol/m 3 [19]. As the degradation and diffusion dominate in the initial time until the monomers are converted into oligomers, the dimensionless autocatalytic concentration reduces exponentially. The concentration profiles have been taken at η = 2 because after that the concentration gets diffused into the extraction medium. The simulation of dimensionless autocatalytic concentration profiles with dimensionless time are shown in Figure 6.  In Figure 7, the dimensionless autocatalytic concentration profiles are presented in the coating film with the distance at different dimensionless time. When τ = 0.1, the dimensionless concentration is at a maximum inside the coating. As the time increases, the concentration in the coating starts to reduce due to the reaction and diffusion phenomenon until = 2 when τ = 0.3. After that, there is an exponential reduction in the concentration as the oligomers are converted to monomers and concentration has dissolved in the extraction medium. This phenomenon can be inferred from the trend at τ = 0.5 and 0.9.
The results for nutrient release profiles and autocatalytic concentration profiles suggest that it is important to consider the reaction-multi diffusion model to study the nutrient release characteristics. The autocatalytic degradation reaction has a significant effect on the release and concentration profiles. This provides further insight to design and develop biodegradable polymer-coated controlled-release fertilizers.  In Figure 7, the dimensionless autocatalytic concentration profiles are presented in the coating film with the distance at different dimensionless time. When τ = 0.1, the dimensionless concentration is at a maximum inside the coating. As the time increases, the concentration in the coating starts to reduce due to the reaction and diffusion phenomenon until = 2 when τ = 0.3. After that, there is an exponential reduction in the concentration as the oligomers are converted to monomers and concentration has dissolved in the extraction medium. This phenomenon can be inferred from the trend at τ = 0.5 and 0.9.
The results for nutrient release profiles and autocatalytic concentration profiles suggest that it is important to consider the reaction-multi diffusion model to study the nutrient release characteristics. The autocatalytic degradation reaction has a significant effect on the release and concentration profiles. This provides further insight to design and develop biodegradable polymer-coated controlled-release fertilizers. In Figure 7, the dimensionless autocatalytic concentration profiles are presented in the coating film with the distance at different dimensionless time. When τ = 0.1, the dimensionless concentration is at a maximum inside the coating. As the time increases, the concentration in the coating starts to reduce due to the reaction and diffusion phenomenon until η = 2 when τ = 0.3. After that, there is an exponential reduction in the concentration as the oligomers are converted to monomers and concentration has dissolved in the extraction medium. This phenomenon can be inferred from the trend at τ = 0.5 and 0.9.

Conclusions
An analytical-numerical solution is presented to the reaction-multi diffusion model for a spherically-shaped PLA-coated CRF granule. We employed the Zakian method for the first time (to the best of the authors' knowledge) for numerical inversion of the Laplace transformation that arises from the reaction-multi diffusion system. The expression gives an intuition of (i) nutrient release profiles and (ii) autocatalytic concentration profiles for different conditions, such as diffusion-only, both reaction and diffusion, and degradation-only reactions for the case of autocatalytic concentration. The limitation of this model occurs when the diffusion coefficient is considered as independent of the autocatalytic concentration profiles and time. Other degradation effects, such as enzymatic and microbial degradation have not been carried out but, in fact, it strongly affects the nutrient release.
The solution indicates that the reaction rate is a key parameter for the prediction of transition between diffusion and degradation of CRF system when autocatalytic concentration and nutrient release profiles are taken into account. The coupling effect of the reaction and diffusion system for nutrient release from PLA-coated controlled-release fertilizer is incorporated in this model. This helps to design and develop biodegradable polymer-coated CRF systems.