CET Volume 86 DOI: 10.3303/CET2186213 Paper Received: 28 August 2020; Revised: 18 March 2021; Accepted: 4 May 2021 Please cite this article as: Bouchkira I., Benjelloun S., Khamar L., Latifi A.M., 2021, Thermodynamic-based Model for the Prediction of the Fouling Phenomena in a Wet Phosphoric Acid Process, Chemical Engineering Transactions, 86, 1273-1278 DOI:10.3303/CET2186213 CHEMICAL ENGINEERING TRANSACTIONS VOL. 86, 2021 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Sauro Pierucci, Jiří Jaromír Klemeš Copyright © 2021, AIDIC Servizi S.r.l. ISBN 978-88-95608-84-6; ISSN 2283-9216 Thermodynamic-Based Model for the Prediction of the Fouling Phenomena in a Wet Phosphoric Acid Process Ilias Bouchkiraa,b,c, Saad Benjelloun b, Lhachmi Khamar b,d, Abderrazak M. Latifia,b,* a Laboratoire Réactions et Génie des Procédés, CNRS-ENSIC, Université de Lorraine, Nancy, France. b Université Mohammed VI polytechnique, Benguerir, Maroc. c LCPCE, Université Hassan II, Casablanca, Maroc d LIPIM, Université Sultan Moulay Slimane, Beni-Mellal, Maroc abderrazak.latifi@univ-lorraine.fr This paper deals with the development of a thermodynamic-based model for the prediction of fouling phenom- ena in a wet phosphoric acid process. The developed model is based on mass and heat balances along with equations of Pitzer's thermodynamic model. It involves several unknown parameters that should be identified from experimental data. For this purpose, a database of experimental measurements carried out in this work, and other measurements from the literature are used. The comparison of the available measurements to the predictions generated using the values of the identified parameters shows the high performance of the devel- oped model. Furthermore, the model is used to predict thermodynamic saturation indices of the main minerals that are likely to precipitate during the digestion process, and may cause fouling. A global sensitivity algorithm is finally used to determine the operating parameters that most influence the precipitation of these mineral phases. The results are very consistent with previous works and with the measurements carried out on indus- trial plants. Moreover, they explain many dysfunctions observed in phosphoric acid production facilities. 1. Introduction The wet phosphoric acid process is one of the widely used processes to produce industrial phosphoric acid. It consists of a series of reactors where a phosphate-containing mineral (usually apatite) is digested using a concentrated sulfuric acid solution. This process is widely used in industry since it is simple to implement and adaptable to various types of phosphate rock. During the digestion reaction, sulfuric acid ( ) dissociates into and ions, and at the same time, the phosphate rock releases many phosphorus and calcium-based ions. These ions contribute to the overall digestion reaction (Eq. (1)) (Becker, 1989) which produces liquid phosphoric acid ( ) and solid calcium sulfates ( ). ( ) + 3 + → 2 + 3 , ; ∈ {0;0.5;2} (1) The reaction shows that ions (coming from sulfuric acid) attack the phosphorus-based ions to produce liq- uid phosphoric acid as the main product, whereas the ions react with the calcium-based ions to produce solid calcium sulfates ( ) as the main by-product. Depending on the operating temperature, the cal- cium sulfates can be either gypsum ( 2 ) , bassanite ( 0.5 ) or anhydrite ( ). The digestion reactor outlet stream is separated into liquid and solid products in a filtration unit located down- stream of the reactors. The liquid phosphoric acid is then recovered for purification and marketing, while the solid calcium sulfate is stored for later treatment and valorization. 1273 Figure 1: Schematic representation of an industrial phosphate digestion tank with three sections and nine re- actors On an industrial scale (e.g., phosphoric acid plants located in Jorf Lasfar, Morocco), the reaction is carried out in a cylindrical tank made-up of several continuous reactors of the same volume and uniformly distributed (Figure 1). The tank is divided into three sections of three reactors each. Sulfuric acid feed flowrate is spitted into three streams each feeding a section. The phosphate ore is fed into the first reactor of the first section and the product leaves the tank at the last reactor of the last section. The side reactions, that occur due to impurities coming from the used reagents, are among the problems en- countered while using the wet phosphoric acid processes. Indeed, these impurities cause the precipitation of mineral deposits leading to fouling phenomena within the reactors, heat exchangers, pumps, etc. These foul- ing problems often lead to production shut down for maintenance which is very costly in terms of time, produc- tivity, and safety. The main impurities can be classified into two categories (Becker, 1989; Azaroual et al, 2012). (i) Impurities favoured by the used reagents such as sulfur, fluorine, silica, magnesium, aluminium, or- ganic elements, and heavy metals, and (ii) impurities favoured by the process itself and consist mainly of so- dium and potassium coming from the process water. Based on literature works and history of the exploitation of the wet phosphoric acid processes, it appears that in addition to the production of gypsum ( 2 ), bassanite ( 0.5 ) and anhydrite ( ) due to the super saturation of the reaction medium by the free sulfates (Becker, 1989), other solid phases can pre- cipitate due to the side reactions caused by the presence of impurities. Indeed, Zmemla et al. (2020), Khamar (2012) and Frazier et al. (1977) proved the existence of brushite ( 2 ) and shukhrovite ( :12 ) by means of scanning electron microscope and X-ray diffraction studies. Furthermore, quantitative X-ray diffraction analysis of quartz are carried out and proved its existence under several crystal- line forms including the cristobalite, chalcedony, ( ) and ( ). They are produced mainly because of the presence of silica, and their crystalline shape depends on the operating conditions, in particular, tempera- ture and pressure (Zmemla et al., 2020; Khamar, 2012; Waerstad, 1985). Moreover, (Zmemla et al., 2020; Skafi et al., 2019; Frayret et al., 2006; Khamar, 2012) analysed the behaviour of malladrite ( ) and hier- atite ( ) using X-ray diffraction techniques and even through thermodynamic modelling based on water activity measurements. These minerals are produced during the digestion due to the use of untreated process water (usually used for dilution, cooling and other purposes) which brings sodium and potassium components that capture ions of the same polarity and opposite charges such as silica and fluorine ions. Finally, magnesi- um-based components, e.g. magnesium hexafluorosilicate ( : , = 0,6,7), are produced during the digestion and may also cause fouling problems. In many works found in literature, these components are studied in order mainly to assess their impact on the filterability of the produced gypsum cake (Rashad et al., 2004; Khamar, 2012). In this work, we are interested in the simultaneous study of all these competitive fouling phenomena. More specifically, we focus on the development of a thermodynamic model to predict the precipitation of the previ- ously described solid phases and to determine their sensitivity with respect to the operating conditions using a global sensitivity analysis approach. 2. Thermodynamic modelling A thermodynamic-based model is developed based on mass and heat balance equations, along with the Pitzer thermodynamic model equations (Pitzer, 2018). It is exploited to predict the saturation indices of solid phases used to determine whether a given mineral can precipitate during the digestion process, thus leading to fouling phenomena. 1274 2.1 Saturation index of minerals The saturation index is used as the main criterion to determine whether the reactive mixture is saturated, un- der-saturated or super-saturated with respect to a given mineral. The saturation index of a mineral that disso- ciates or associates as … ⇌ + ⋯+ is calculated as: = log = ( ) (2) where is the ion activity product, is the solubility product of the mineral, is the unknown molality of component and is its activity coefficient, is the stochiometric coefficient. 2.2 Thermodynamic-based model The developed thermodynamic-based model aims to compute the saturation indices of the considered miner- als (listed in Table 1) given the operating conditions. As shown in Eq. (2), to compute the saturation indices of these minerals, the molalities and the activity coefficients of 25 aqueous elements (Table 1) are need- ed. They are deduced from a set of equations consisting of 9 material balances (Eqs. (4)), 1 charge balance (Eq. 5), 15 equilibrium constants (presented in Table 2) calculated using Eqs. 6 and 7, and 23 activity coeffi- cients calculated using the Pitzer thermodynamic model (Pitzer, 2018). Finally, it is noteworthy that a heat bal- ance is supported by the model to calculate the operating temperature (Eq. 8). The resulting nonlinear system of equations involves 50 unknown variables, i.e., 25 molalities and 25 activity coefficients, and 50 equations, which can be solved using a suitable numerical solver. ( ) = , ∈ { , , , , , , , , } (4) = 0 (5) = ( ) , (6) log = log + ∆ 1 − 1 , = 1,…, (7) + + = + (8) where and are the number of components and reactions respectively, ( ) is the given total inlet concentration of element in Eq. (1), is the molality of component involved in the reactions listed in Ta- ble 1, , is equal to the number of atoms of element in component . For example, the mass balance on (F) develops as: ( ) = 6 H2SiF6 + F− + 6 SiF6 + 4 AlF4 + 6 . are the electrical charges, is the activity coefficient of component . , is the stoichiometric coefficient of component involved in reaction . and ∆ refer to the equilibrium constant and the enthalpy of the reaction at = 298 respectively, their values are reported in Table 2. , , , and are the enthalpy of the inlet products, the agitation heat, the sulfuric acid dilution heat, the enthalpy of the slurry leaving the reactor and the heat to be removed from the process to control the operating temperature respectively. The Pitzer thermodynamic model (Pitzer, 2018) is used to compute the activity coefficient as: ln( ) = 2 + 2 , + ,, + 3 , ,, (9) where i, j,k refer to different components and is the Debye-Huckel function. λ , , λ , and ψ , , are the un- known Pitzer model parameters, they should be determined from the equilibrium measurements carried out in this work and presented in the next section. Furthermore, a suitable estimability analysis method should be used to determine the most estimable parameters from the available data. 1275 Table 1: Mineral phases and their formation-dissolution reactions supported by the developed model Mineral Reactions (S ⇌ Aq) log( ) References Fluorapatite Ca (PO ) F + H ⇌ HF + 3 PO + 5 Ca -0.910 (Gustafsson, 2011; Ball and Nordstrom, 1991) Gypsum CaSO :2H O ⇌ Ca + SO + 2 H O -4.848 (Gustafsson, 2011) Bassanite CaSO :0.5H O ⇌ Ca + SO + 0.5 H O -3.661 (Gustafsson, 2011) Anhydrite CaSO ⇌ Ca + SO -4.637 (Gustafsson, 2011) Brushite CaHPO :0.5H O ⇌ Ca + HPO + 0.5H O -6.550 (Ball and Nordstrom, 1991) Shukhrovite Ca SO AlSiF ⇌ CaSO + SiF + AlF + 3CaF -19.55 (Zmemla et al., 2020; Fraizer, 1977) Malladrite Na SiF ⇌ SiF + 2 Na -6.475 (Skafi et al., 2019 ; Azaroual et al. 2012) Hieratite K SiF ⇌ SiF + 2 K -7.039 (Skafi et al., 2019) Magadiite NaSi O (OH) :3H O + H + 9H O ⇌ Na + 7H SiO -14.30 (Gustafsson, 2011; Ball and Nordstrom, 1991) Fluorite CaF ⇌ Ca + 2 F -10.96 (Gustafsson, 2011; Ball and Nordstrom, 1991) Quartz SiO + 2 H O ⇌ H SiO -3.018 (Gustafsson, 2011; Ball and Nordstrom, 1991) MgSiF6 MgSiF ⇌ SiF + Mg -9.960 (Ball and Nordstrom, 1991; Khamar, 2012) MgSiF6.6H2O MgSiF :6H O ⇌ SiF + Mg + 6 H O -12.96 (Ball and Nordstrom, 1991; Khamar, 2012) MgSO4.7H2O MgSiF :7H O ⇌ SiF + Mg + 7 H O -14.96 (Ball and Nordstrom, 1991; Khamar, 2012) CaSiF6 CaSiF ⇌ Ca + SiF -0.519 (Gustafsson, 2011; Ball and Nordstrom, 1991) Table 2: Aqueous phases supported by the developed model Equilibria Reactions ∆ (J/mol) References 1 HSO ⇌ SO + H 0.0103 -1.274. 10 (Gustafsson, 2011; Pitzer, 2018) 2 H PO ⇌ PO + H 0.0071 -4.271. 10 (Gustafsson, 2011; Ball and Nordstrom, 1991) 3 PO ⇌ HPO + H 6.3. 10 -1.800. 10 (Gustafsson, 2011; Ball and Nordstrom, 1991) 4 HPO ⇌ PO + H 4.2. 10 -1.500. 10 (Gustafsson, 2011; Ball and Nordstrom, 1991) 5 H PO + PO ⇌ 0.2550 -9.266. 10 (Khamar, 2012) 6 HF ⇌ F + H 7.2.10 -1.329. 10 (Gustafsson, 2011) 7 H SiF ⇌ SiF + 2H 0.3. 10 -4.799. 10 (Gustafsson, 2011) 8 H SiO ⇌ SiO + H O 3.52.10 -1.121. 10 (Gustafsson, 2011) 9 CaF ⇌ Ca + F 1.148. 10 -1.255. 10 (Gustafsson, 2011; Ball and Nordstrom, 1991) 10 CaSO ⇌ Ca + SO 0.6394 -7.200. 10 (Pitzer, 2018; Shen et al., 2020) 11 AlSO ⇌ Al + SO 9.549. 10 -8.368. 10 (Gustafsson, 2011; Ball and Nordstrom, 1991) 12 AlF ⇌ Al + 4F 9.772. 10 -8.368. 10 (Gustafsson, 2011; Ball and Nordstrom, 1991) 13 MgSO ⇌ Mg + SO 5.623. 10 -4.184. 10 (Gustafsson, 2011; Ball and Nordstrom, 1991) 14 NaSO ⇌ Na + SO 0.1995 -4.170. 10 (Gustafsson, 2011; Ball and Nordstrom, 1991) 15 KSO ⇌ K + SO 1.1412 -8.368. 10 (Gustafsson, 2011; Ball and Nordstrom, 1991) 3. Experimental measurements Several equilibrium measurements of sulfuric and phosphoric acid solutions are carried out in the laboratory, they consist of pH, conductivity, and density and were used to calculate the molalities of the different ions and molecules (i.e., speciation) involved in the reactive mixtures. To perform these measurements, 98w% sulfuric acid and 65w% phosphoric acid are used. Sulfuric acid measurements are carried out for concentrations rang- ing from infinite dilution to 4 moles/kgw and temperatures ranging from 298 K to 353 K, whereas, measure- ments for phosphoric acid are performed for concentrations varying between infinite dilution and 12 moles/kgw, and in the same range of temperature. Several experimental data are also collected from the liter- ature (Shen et al., 2020), they consist of the solubilities of gypsum and anhydrite calcium sulfates in sulfuric acid, carried out at different molality and temperature conditions. 4. Results and discussion Experimental measurements are used to identify the unknown parameters of the Pitzer thermodynamic model. Beforehand, a parameter estimability analysis was carried out based on our recently developed global estima- bility algorithm (Bouchkira et al., 2021). This algorithm is used to determine the unknown parameters that are estimable from the available measurements. The most estimable parameters are then identified by minimizing the square-error between the model predictions and the available experimental measurements. Furthermore, the optimization problem was implemented within GAMS environment and a branch-and-reduce algorithm was used to reach the global minimum. The computation time was about three hours using a Dell Precision T7810 Bi-Xeon 12 x Core 64GB Dell Precision workstation. The non-estimable parameters and those for which no experimental measurements are available for their identification are fixed from the literature and from previous works. The values of the interaction parameters are used to compare the predictions of the developed model with experimental measurements, and thus to evaluate the model accuracy. Figures 2 A-C present an exam- ple of results, they show the prediction of the speciation of sulfuric and phosphoric acid solutions as well as the solubility of gypsum and anhydrite calcium sulfates at 298 K. The Pearson product-moment correlation co- efficient is also calculated and reported for each case. Its high values demonstrate the very good perfor- mance of the developed model. 1276 Figure 2: (A) and (B): Sulfuric and Phosphoric acids speciation, (C): gypsum and Anhydrite solubility Figure 3: (A): Saturation indices of minerals, (B): Total sensitivity with respect to operating conditions The saturation indices of minerals in Table 1 are calculated using the developed model and shown in Fig. 3A. Moreover, a global sensitivity analysis algorithm developed by Saltelli (2002) is used to assess the sensitivity of each saturation index with respect to the most important operating conditions (i.e., sulfur and phosphorus concentrations and temperature (Fig. 3.B)). The results show that fluorapatite, hydroxyapatite, fluorite and cal- cium hexafluorosilicate have negative saturation indices meaning that the thermodynamic conditions of the digestion are favorable for their dissolution and not for their crystallization. This is consistent with the results of Becker (1989) and who studied the dissolution of phosphate-containing minerals in phosphoric acid production processes. The saturation indices of the calcium-based minerals (i.e., gypsum, bassanite and anhydrite) are positive meaning that these minerals are likely to crystallize. Indeed, during the digestion process, the calcium ions released from the phosphate ore tend to capture the free sulfate ions of the same polarity and opposite charges leading to calcium sulfates crystallization. The sensitivity analysis results show that the saturation in- dices of these minerals are very sensitive to the operating temperature which is very consistent with previous studies (Becker, 1989; Shen et al., 2020). In general, these minerals are considered as principal products of the digestion process and are not involved in the fouling of the process facilities. On the other hand, it is shown that the operating conditions of the digestion are favorable to produce the brushite and shukhrovite as proved by Fraizer (1977) and Zmemla et al., (2020) who studied the existence of these minerals in the gyp- sum cake obtained in the filtration unit. Meanwhile, Khamar (2012) proved the crystallization of these mineral phases during the cooling of industrial phosphoric acid in storage tanks. Sensitivity analysis shows that the production of these minerals depends mainly on sulfur concentration which is consistent with the results of Becker (1989) and Zmemla et al., (2020). The developed model predicts positive saturation indices for mal- ladrite, hieratite and magadiite, which are mainly produced from the combination of sodium and potassium fa- voured using untreated process water, and silica and fluorine from the phosphate ore. Azaroual et al. (2012) showed in their work that malladrite presents 68% of the products that cause fouling of the concentration units downstream of the digestion tank, and 96% of the products that cause fouling of the heat exchangers. It is al- so shown that quartz tends to precipitate in different crystalline forms like chalcedony, cristobalite, ( ), ( ). Recently Khamar (2012) and Zmemla et al. (2020) showed their existence by means of X-ray and scanning electron microscopy technologies and their high sensitivity to the operating temperature, moreover, to limit their precipitation, the industrialists advise doping the reaction mixture with aluminum, which favors the fixation of silica and the production of clay and shukhrovite, which are favorable for improving the filterability of the gypsum cake in view of their crystalline form. Finally, the saturation indices of magnesium-containing min- 1277 erals are positive and depend mainly on sulfur concentration. In fact, magnesium is one of the critical impuri- ties involved in several side reactions which cause fouling phenomena and decrease the filterability of the gypsum cake. It is noteworthy that these magnesium-containing minerals are generally characterized by very small sizes, which favours the fouling phenomena (Khamar, 2012) which minimizes the porosity of the gypsum cake and degrades its filterability (Becker, 1989). 5. Conclusions In this work, a thermodynamic model is developed to enable the prediction of minerals that may precipitate during the digestion of phosphate ore due to side reactions and lead to fouling phenomena. It has been shown that under current production conditions, malladrite, hieratite, magadiite, quartz and magnesium hexafluorosil- icate in their different crystalline forms are the most critical minerals that are likely to precipitate and cause fouling phenomena. It is also shown that, limiting the precipitation of some of them requires mainly to operate under an optimal sulfur concentration and temperature which allows to neutralize them or even to limit their production. Moreover, the treatment of the used process water would be an effective solution that should be implemented by the phosphate industries to minimize these fouling problems. Other side reactions that take place during the digestion process were neglected for the sake of simplification; they are mainly due to heavy metals, organic compounds, etc. They should be considered to further improve the accuracy of model. References Azaroual, M., Kervevan, C., Lassin, A., André, L., Amalhay, M., Khamar, L., El Guendouzi, M., 2012, Thermo- kinetic and physico-chemical modeling of processes generating scaling problems in phosphoric acid and fertilizers production industries, Procedia Engineering, 46, 68-75. Ball, J. W., Nordstrom, D. K., 1991, WATEQ4F--User's manual with revised thermodynamic data base and test cases for calculating speciation of major, trace and redox elements in natural waters (No. 90-129). Becker, P., 1983, Phosphates and phosphoric acid. Marcel Dekker, Inc. Bouchkira, I., Latifi, A. M., Khamar, L., Benjelloun, S., 2021, Global sensitivity based estimability analysis for the parameter identification of Pitzer’s thermodynamic model, Reliability Engineering & System Safe- ty, 207, 107263. Cherif, M., Mgaidi, A., Ammar, M. N., Abderrabba, M., Fürst, W., 2002, Representation of VLE and liquid phase composition with an electrolyte model: application to H3PO4–H2O and H2SO4–H2O, Fluid phase equilibria, 194, 729-738. Frayret, J., Castetbon, A., Trouve, G., Potin-Gautier, M., 2006, Solubility of (NH4)2SiF6, K2SiF6 and Na2SiF6 in acidic solutions, Chemical physics letters, 427(4-6), 356-364. Frazier, A. W., Lehr, J. R., Dillard, E. F., 1977, Chemical behaviour of fluorine in production of wet-process phosphoric acid, Environmental Science & Technology, 11(10), 1007-1014. Gustafsson, J. P., 2011, Visual MINTEQ 3.0 user guide. KTH, Department of Land and Water Recources, Stockholm, Sweden. Khamar, L., 2012, Experimental study and thermodynamic modeling of fouling problems in industrial phos- phoric acid production facilities. Thesis report. University of Hassan II. Casablanca. Pitzer, K. S., 2018, Activity coefficients in electrolyte solutions, CRC press. Saltelli, A., 2002, Making best use of model evaluations to compute sensitivity indices, Computer physics communications, 145(2), 280-297. Rashad, M. M., Mahmoud, M. H. H., Ibrahim, I. A., Abdel-Aal, E. A., 2004, Crystallization of calcium sulfate dihydrate under simulated conditions of phosphoric acid production in the presence of aluminum and mag- nesium ions, Journal of crystal growth, 267(1-2), 372-379. Shen, L., Sippola, H., Li, X., Lindberg, D., Taskinen, P., 2020, Thermodynamic Modeling of Calcium Sulfate Hydrates in a CaSO4–H2SO4–H2O System from 273.15 to 473.15 K up to 5 m Sulfuric Acid, Journal of Chemical & Engineering Data, 65(5), 2310-2324. Skafi, M., Elyamani, Y., & El Guendouzi, M., 2019, Solubility of Mixed Ternary Systems of Hexafluoridosilicate Salts (M= Na+, K+, or NH4+) in Hexafluoridosilicic Acid Aqueous Solutions at T= 353.15 K. Journal of Chemical & Engineering Data, 64(8), 3290-3299. Waerstad, K. R., 1985, Quantitative X-ray Diffraction Analysis of Calcium Sulfates and Quartz in Wet-Process Phosphoric Acid Filter Cakes. Advances in X-ray Analysis, 29, 211-216. Zmemla, R., Sdiri, A., Naifar, I., Benjdidia, M., Elleuch, B., 2020. Tunisian phosphogypsum tailings: Assess- ment of leaching behavior for an integrated management approach, Environmental Engineering Re- search, 25(3), 345-355. 1278