Computers and Chemical Engineering 156 (2022) 107536

Contents lists available at ScienceDirect

Computers and Chemical Engineering

journal homepage: [Link]/locate/compchemeng

Modeling and multi-objective optimization of the digestion tank of an

industrial process for manufacturing phosphoric acid by wet process
Ilias Bouchkira a,b,c, Abderrazak M. Latifi a,b,∗, Lhachmi Khamar a,d, Saad Benjelloun b
a
Laboratoire Réactions et Génie des Procédés, CNRS-ENSIC, Université de Lorraine, Nancy Cedex, France
b
Mohammed VI Polytechnic University, Lot 990 Hay Moulay Rachid, Benguerir, Morocco
c
Laboratory of Physical-Chemistry, Catalysis & Environment, University Hassan II-Casablanca, Casablanca, Morocco
d
Laboratoire d’Ingénierie des Procédés, Informatique et Mathématiques, Université Sultan Moulay Slimane, Beni-Mellal, Morocco

a r t i c l e i n f o a b s t r a c t

Article history: This paper deals with the modeling and multi-objective optimization of an industrial phosphoric acid
Received 28 February 2021 process. The objective is to determine the operating conditions which minimize the chemical losses of
Revised 6 August 2021
phosphate and maximize the productivity of the digestion tank. To achieve these objectives, a process
Accepted 9 September 2021
model, based on mass, charge, and energy balances along with thermodynamic equilibrium equations is
Available online 12 September 2021
developed. Experimental measurements concerning sulfur and phosphorus based systems are carried out
Keywords: and other measurements concerning fluorine, silicon and calcium sulfates based systems are collected
Phosphoric acid process from the literature. The results show that the predictions are in very good agreement with the mea-
Digestion tank surements. The developed model is then used in a multi-objective optimization problem of an industrial
Modeling manufacturing process to determine the set of optimal operating conditions. The optimization problem
Multi-objective optimization is solved by means of epsilon-constraint method and the optimal solutions are ranked using the multi-
Surrogate model
attribute utility theory. The best solutions are compared to industrial measurements based on surrogate
Decision making
modeling and are found to be consistent with the current operating conditions. Their implementation
would significantly improve the current process performances.
© 2021 Elsevier Ltd. All rights reserved.

1. Introduction each feeding one of the three sections. The phosphate ore is fed
into the first reactor of the first section and the final product leaves
Phosphoric acid is one of the most produced and marketed the tank at the last reactor of the third section (Fig. 1).
chemicals in the world. It is mainly used in the manufacture of fer- During the digestion process, as shown by the overall reaction
tilizers (Rickard, 20 0 0; Li et al., 2016; Leikam and Achorn, 2005), below, sulfuric acid is dissociated into H + and SO24− ions. The H +
in the food industry (Amin et al., 2010; Lampila, 2013) and even in ions are involved in the extraction of the phosphate elements and
the pharmaceutical industry (Geeson and Cummins, 2018). In sev- their diffusion from the solid phase (phosphate ore) towards the
eral producing countries, it is produced by the wet process which liquid phase (phosphoric acid) whereas the SO24− ions crystallize
consists of reacting the phosphate ore with sulfuric acid (Becker, with the calcium ions Ca2+ to form a solid calcium sulfate which
1989; Slack and James, 1973; Dorozhkin, 1996). can be either anhydrite C aSO4 , gypsum C aSO4 :2H2 O or bassanite
On industrial scale, phosphoric acid is mainly produced by the CaSO4 :0.5H2 O depending on the operating conditions, in particu-
digestion of the phosphate ore using a concentrated sulfuric acid lar, the temperature and the sulfates concentration (Becker, 1989;
solution. The digestion is carried out in a cylindrical tank which Slack and James, 1973; Dorozhkin, 1996). The liquid phosphoric
consists of a series of nine continuous reactors of the same vol- acid and solid calcium sulfates are separated in a vacuum filtration
ume and uniformly distributed inside the tank. The latter is further unit downstream. Phosphoric acid is recovered for valorization and
divided into three sections of three continuous reactors each. The marketing, while solid calcium sulfates are eliminated as an unde-
sulfuric acid feed flow rate is therefore split into three streams, sirable product.
Ca3 (PO4 )2 + 3H2 SO4 + 3x. H2 O −→ 2H3 PO4 + 3CaSO4 , x H2 O

∗
Moreover, the production of more gypsum enhances the produc-
Corresponding author at: Laboratoire Réactions et Génie des Procédés, CNRS-
ENSIC, Université de Lorraine, Nancy Cedex, France. tivity by increasing the phosphoric acid production, whereas the
E-mail addresses: [Link]@[Link] (I. Bouchkira), production of the bassanite increases the viscosity of the reactive
[Link]fi@[Link] (A.M. Latifi). mixture within the reactors and consequently lowers the perfor-

[Link]
0098-1354/© 2021 Elsevier Ltd. All rights reserved.
I. Bouchkira, A.M. Latifi, L. Khamar et al. Computers and Chemical Engineering 156 (2022) 107536

Fig. 1. Schematic representation of an industrial wet phosphoric acid process.

mances of the filtration process downstream. It is therefore very The dihydrate processes (studied in this work) are very popular
important to control the operating conditions to favor the produc- and widely used in industry. They are simple to implement, can
tion of gypsum and limit the formation of the bassanite. be adapted to various types of minerals, and have the advantage
On the other hand, the required amount of sulfuric acid should to be easily controllable (Gilmour, 2013). The produced solid phase
be optimally used in order to optimize the production perfor- in these processes consists mainly of gypsum. Regarding the hemi-
mance. Indeed, a deficit in sulfuric acid during the digestion causes hydrate processes, they are designed to produce a more concen-
a decrease in sulfate ions SO24− concentration in the digestion tank. trated phosphoric acid, which is advantageous in terms of energy
Therefore, the free calcium ions (coming from the raw rock) tend consumption since the existence of the concentration unit down-
to capture the phosphate ions HP O24− and crystallize them as solid stream of the digestion, which is energy intensive, is no longer
brushite CaHP O4 :2H2 O. The latter decreases the chemical yield of justified. The main disadvantages of the hemihydrate processes
the process since the solid phase is discharged as undesirable are due to the higher reaction temperature, which causes corro-
product leading to phosphate losses by syncrystallization (Becker, sion and foaming problems. Moreover, the produced solid phase in
1989; Dorozhkin, 1996). these processes consists of bassanite which lowers the filterability
On the contrary, the use of an excess of sulfuric acid causes not performances and leads to some phosphate losses during the fil-
only a local increase of the temperature due to the heat of dilu- tration (Cho et al., 1996). To ovoid these disadvantages, the hemi-
tion, but also leads to the formation of a gypsum layer around the hydrate with re-crystallization processes are designed. They proved
phosphate particles (particle coating) thus protecting them from to be efficient in terms of the concentration of the produced phos-
further attack by the sulfuric acid and resulting also in phosphate phoric acid and the filterability of the calcium sulfate produced.
losses (unattacked losses). The appropriate amount of sulfuric acid However, their control still poses many problems compared to the
along with its optimal distribution over the three sections of the dihydrate and hemihydrate processes (Theys et al., 2014). The com-
digestion tank are then relevant to the optimal operation of the bination of the dihydrate and hemihydrate processes are also used
digestion tank. (Alfariss et al., 1992; Gurgul et al., 2019), they are based on the
In this paper, the objective is to determine the operating con- production of bassanite and then its transformation into gypsum.
ditions of the digestion tank that minimize the phosphate ore The advantage of this combination of the two technologies is that
losses (syncrystallized and unattacked losses) while meeting the they are less expensive, since their design and implementation re-
constraints on the excess of sulfuric acid above the stoechiomet- quire a lower capital investment and especially lower operation
ric amount and on the reaction temperature. Beforehand, Pitzer’s cost. However, their disadvantages are mainly due to the corrosion
thermodynamic model, needed in optimization, is calibrated using problems and the complexity of their control.
experimental measurements. The novelty of the paper is threefold: Finally, dehydrate to hemihydrate processes have very similar
(i) global estimability analysis and identification of the most es- advantages and disadvantages, the only difference is that the reac-
timable parameters of Pitzer’s thermodynamic model from experi- tion temperature is lower, which avoids the corrosion and foam-
mental data, (ii) exploitation of this model for the multi-objective ing problems. However, these processes are designed to produce
optimization of an industrial phosphoric acid manufacturing pro- gypsum and its transformation into bassanite, which reduces the
cess and selection of the best optimal solution using a decision filtration performances.
making-aid method, (iii) and experimental validation on an indus- Modeling and simulation of the wet phosphoric acid process
trial scale of the optimal solution. have been addressed in the literature for both hemihydrate (Al-
Fariss et al., 1993; Cho et al., 1996; Gioia et al., 1977; Grema
2. Literature review on wet process et al., 2018; Yang et al., 2018) and dehydrate processes (Al-Thyabat
and Zhang, 2015; Bharathi et al., 2013; Peng et al., 2015; Rabad-
In industry, there is a wide variety of "wet processes" designed jieva et al., 2020; Bichri et al., 2020). However, their optimal de-
to produce liquid phosphoric acid and solid calcium sulfates un- sign and operation have not been sufficiently investigated and
der specific operating conditions. They mainly include : (i ) dihy- very few research contributions are devoted to the optimization
drate, (ii ) hemihydrate, (iii ) hemihydrate to dihydrate, (iv ) hemi- issues either of the entire process or of one or more of its
hydrate with re-crystallization, (v ) and dihydrate to hemihydrate unit operations. Among these works, Gioia et al. (1977) consid-
processes (Abu-Eishah and Abu-Jabal, 2001; Theys et al., 2014; Bel- ered the entire hemihydrate process and determined the struc-
boom et al., 2015). ture of the process that minimizes an economical performance

2
I. Bouchkira, A.M. Latifi, L. Khamar et al. Computers and Chemical Engineering 156 (2022) 107536

index which is the sum of annual fixed and operating costs.

1
 j
Papadopoulos and Seferlis (2009) also optimized a complete indus-
f1 = log aCa2+ .aHPO2− .a2H2 O .k−1
s1 (1)
trial process for the production of phosphoric acid by determining 3 4
j
the configuration that maximizes the amount of acid produced.
Later, Joao et al. (2016) focused on a phosphoric acid concentra- 1
 j
tion unit and used dynamic simulation tools to study the opti- f2 = log aCa2+ .aSO2− .a2H2 O .k−1
s2 (2)
3 4
mization of its operability. Recently, Ma et al. (2018) improved the j

flowsheet of the wet-process phosphoric acid production to mini- where ai is the activity of component i in the reaction mixture,
mize the flow rate of the fresh water introduced into the process. ks1 and ks2 are the solubility products of the brushite and of the
Bouchkira et al. (2021a) developed a thermodynamic-based model gypsum respectively and, j refers to the number of each section.
for the prediction of fouling phenomena in the wet phosphoric
acid manufacturing process. Multi-objective optimization was used
3.2. Decision variables
by Hemalatha and Rani (2017) to develop optimal operating recipes
in a batch cooling reactor for unseeded and seeded crystallization
The decision variables consist of the excess of sulfuric acid Sc
processes. Very recently, (Bouchkira et al., 2021b; 2021d) carried
over the stoichiometric amount required by the reactions, its distri-
out a multi-objective optimization of the digestion tank of an in-
bution ratios over the sections, i.e., w1 , w2 and w3 , and the cooling
dustrial phosphoric acid manufacturing process.
heat Qcool to be evacuated to control the temperature of the reac-
tion (Fig. 1).
3. Optimization problem formulation
3.3. Equality and inequality constraints
A multi-objective optimization problem is formulated and
solved in this work in order to determine the trade-offs between Two inequality constraints are introduced to account for the
the following objectives: (i) limit the syncrystallized losses by min- last two objectives. The first one sets the lower and upper lim-
imizing the precipitation of the solid brushite, (ii) enhance the pro- its of sulfuric acid excess to specific values which minimize
ductivity of the process by maximizing the production of gypsum, the unattacked losses. These limits correspond to 1% and 3%
(iii) limit the unattacked losses by setting the excess of sulfuric (Becker, 1989). The second inequality constraint addresses the op-
acid over the stoichiometric amount between two limit values that erating temperature that should remain below 355 K to avoid the
ensure the optimal operation of the process, (iv) control the tem- production of the solid bassanite. These inequality constraints are
perature of the reaction in order to limit the production of the bas- expressed as :
sanite and thus improve the filtration downstream of the process.
1 % ≤ Sc ≤ 3 % (3)
The first two objectives are taken into account by means of two
criteria, whereas the remaining two objectives are considered as
constraints. Indeed, the first two objectives can be predicted di- T ≤ 355 K (4)
rectly through thermodynamic modeling, which is not the case for
the last two objectives. For example, the thermodynamic model The equality constraints consist of process model equations which
developed cannot be used to predict the amount of particles that are based on mass balance Eq. (5)), charge balance (Eq. (6)), equi-
will undergo the coating phenomenon which will protect them librium constants equations (Eqs. (7-(8)), and finally heat balance
from further attack of sulfuric acid (objective iii). The model also (Eq. (10)) (The detailed model equations are given in the suple-
cannot be used to predict the filterability of the produced gypsum mentary material at the end of the paper).
cake (objective iv). However, previous experiences in the literature 
NC
have allowed us to consider these two goals as constraints. Fur- (M )tot = δM,i mi , M ∈ {S, P, F , Si, Ca} (5)
thermore, the experimental work of Becker (1989) suggests that i=1
the optimal excess of sulfuric acid lies between 1% and 3%. He
proved that by operating above 3% excess sulfuric acid, the sec- 
NC
zi .mi = 0 (6)
ondary nucleation of the gypsum on the surface of the phosphate
i=1
particles is very favorable and leads to the formation of a coating
layer around them, which protects them from sulfuric acid. More- 
NC
α 
NC

over, numerous reviews of the literature have shown that an excess Kj = ai i j = (mi .γi )αi j , j = 1, ., NR (7)
of 1% sulfuric acid can be considered as the minimum amount to i=1 i=1

H j  1 
cause supersaturation of the reaction medium and then to precip-
1
itate the gypsum (Becker, 1989; Slack and James, 1973; Dorozhkin, log(K j ) = log(K j0 ) + − (8)
1996). Finally, a last inequality is used to account for the last objec- R T T0
tive, and suggests to operate below 355 K to avoid crystallization of where NC is the number of components, NR corresponds to the
bassanite, which is characterized by high viscosity and small size number of reactions, (M )tot is the total concentration of M. mi
particles, thus decreasing the filterability of the gypsum cake. The is the molality of each component i involved in the equilibria in
ingredients of the multi-objective optimization problem considered Table 1, δM,i is the number of element M in component i. For ex-
are then presented in the next section. ample, for phosphorus element δP,H P O− = 2. As an example, the
5 2 8
mass balance on (P ) develops as:
3.1. Objective functions (P )tot = mH3 PO4 + mH2 PO−4 + mHPO42− + 2.mH5 P2 O−8 (9)

The two (conflicting) optimization criteria are defined by the zi are the electrical charges, γi and ai are the activity coefficient
saturation index (Chidambaram et al., 2011) of the brushite to be and the activity of component i respectively. αi j is the stoichiomet-
minimized all over the sections of the process to limit its produc- ric coefficient of component i involved in reaction j. K j0 and H j
tion, and the saturation index of gypsum to be maximized to en- refer to the equilibrium constant and the enthalpy of the reaction
hance the process productivity. These two criteria are expressed as: j at T0 = 298 K respectively, their values are reported in Table 1.

3
I. Bouchkira, A.M. Latifi, L. Khamar et al. Computers and Chemical Engineering 156 (2022) 107536

Table 1
Reactions involved in the digestion tank and their corresponding equilibrium constant and enthalpy at
298K. (-) means that the reaction is total and does not involve an equilibrium.

Equilibria K j0 H j (J/moles ) Database

1 H2 SO4 −→ HSO−
4 +H
+
- - -
K1
2 HSO− 2−
4 ←→ SO4 + H
+
0.0103 -1.2.104 Gustafsson (2011);Pitzer (2018)
K2
3 H3 PO4 ←→ H2 PO−
4 +H
+
0.0071 -4.2.103 Gustafsson (2011)
K3
4 H2 PO− 2−
4 ←→ HPO4 + H
+
4.2.10−12 -1.5.104 Gustafsson (2011)
K4
5 H2 PO− −
4 + H3 PO4 ←→ H5 P2 O8 0.2550 -9.2.103 Gustafsson (2011)
K5
6 HF ←→ F− + H+ 7.2.10−4 -1.3.104 Ball and Nordstrom (1991)
K6
7 HF + F− ←→ HF−
2 5.5000 -1.7.104 Ball and Nordstrom (1991)
K7 2−
8 H2 SiF6 ←→ SiF + 2 H+ 0.3.102 -6.7.104 Gustafsson (2011)
K8
9 Ca2+ + SO24− ←→ CaSO4 0.6394 -7.2.103 Gustafsson (2011), Pitzer (2018)

Finally, the heat balance for each section of the process is ex- Its use depends mainly on the chemical species present in solution.
pressed as: The general expression of this model provides the excess Gibbs en-
ergy for a solution containing nw kg of solvent through the follow-
QH + QAG + QR+D = Qcool + Qout + Qloss (10)
ing relation:
where QH , QAG , QR+D , Qout , Qloss and Qcool refer respectively to  
Gex
the enthalpy of the inlet reactants, the agitation heat, the sulfuric = nw f + mi m j λi, j + mi m j ψi, j,k + . . . (16)
Wn RT
acid dilution heat, the enthalpy of the slurry leaving the reactor, i, j i, j,k
the heat losses of the digestion tank and the heat to be removed
λi, j represents the short-distance binary interactions between
to control the temperature. The enthalpy of the inlet and outlet
solute species i and j. ψi, j,k represents the ternary interaction pa-
streams are expressed as follows:
rameters. R is the perfect gas constant, Wn is the mass of the
QH = Qsa + Q ps + Qr ph (11) solvant, T is the temperature, f is the Debye-Huckel function
(Pitzer, 2018) which depends on the ionic strength I as:

Qout = Vout ρout C pout (Tout − T0 ) (12) 4Aφ I  

f =− ln 1 + I0.5 (17)
b
where:
1
Qsa = msaC psa (Tsa − T0 ) (2π N0 dw /10 0 0 ) 2 (e2 DkT ) 2
1 3
(13) Aφ = (18)
3
where Aφ is the Debye-Huckel constant, b is a characteristic con-
Q ps = Vps ρ psC p ps (Tps − T0 ) (14)
stant. N0 is the Avogadro number, dw is the density of water, D is
the dielectric constant of water at temperature T , k is Boltzmann’s
Qr ph = Vr ph ρr phC pr ph (Tr ph − T0 ) (15) constant and e is the electrical charge. The activity coefficient γi is
given by the Pitzer model as (Pitzer, 2018):
msa , C psa , Tsa are the mass flow rate, the heat capacity and the
temperature of sulfuric acid respectively. Vps , C p ps , ρ ps , Tp are the zi2   
ln(γi ) = f +2 m j λi, j (I ) + zi2 m j mk λ j,k (I )
volume flow rate, the heat capacity, the density and the tempera- 2
j j
ture of the phosphate slurry respectively. Vr ph , C pr ph , ρ ph , Tr ph are 
the volume flow rate, the heat capacity, the density and the tem- +3 mi m j mk ψi, j,k + . . . (19)
perature of the recycled phosphoric acid respectively. Vout , C pout , j,k
ρout , Tout are the volume flow rate, the heat capacity, the density where:
and the temperature of the outlet flow respectively. T0 is the refer-
1 dλi, j ni
ence temperature. It is worth mentioning that the correlations de- I= mi zi2 , λi, j = , mi = (20)
veloped by Becker (1989) are used to estimate the agitation heat 2 dI nw
i
and the heat released by the dilution of sulfuric acid. It should be
In Eq. (19), it is necessary to use some trade-offs when de-
noted that Qloss is greater for small units than for large ones due
ciding about the number of terms to account for. As suggested
to their large surface/volume ratio. Its value is often too low to be
by Pitzer (2018), this number should be chosen depending on the
taken into account and is therefore neglected.
value of the ionic strength of the solutions. Thus, when it comes
Solving the formulated multi-objective optimization problem
to thermodynamic systems with very low ionic strength, only the
first requires the selection and use of an appropriate electrolytic
first term (i.e., the term of Debye Huckel) can be used. When the
model to predict several necessary thermodynamic properties,
ionic strength is not very high, a very good prediction of the nec-
among which, the activities of the elements considered to predict
essary thermodynamic data can be obtained with only binary in-
the values of the objective functions in Eqs. (1-2), as well as the
teraction parameters. Finally, for very concentrated solutions (with
activity coefficients to predict the equilibrium constants in Eq. (7).
high ionic strength), which is the case in this work, binary and
For this purpose, the Pitzer model is selected and used as one of
ternary interaction parameters must be taken into account.
the most accurate and studied models in the literature. This model
Pitzer gives more appropriate forms of his theory by using βi, j
is presented and detailed in the following section. 
to express the parameters λi, j and λi, j . Moreover, these parameters
3.4. Pitzer’s thermodynamic model are temperature dependent (Simoes et al., 2017) and are usually
expressed as Eqs. (21)–(25).
Pitzer’s thermodynamic model (Pitzer, 2018) is widely used to 2
 √ √
describe the interactions between dissolved ions and the solvent. λi, j (I ) = βi,(0j ) + βi,(1j ) 1 − ( 1 + α I ) e −α I (21)
b

4
I. Bouchkira, A.M. Latifi, L. Khamar et al. Computers and Chemical Engineering 156 (2022) 107536

Table 2
Recent experimental measurements from the literature for the identification of the unknown model parameters and
the validation of the models results. Molality is expressed in (moles/kgw ) and temperature in (K ).

System Type of data Molality Temperature Reference

H3 PO4 − H2 O Osmotic coefficient 0–10.0 298–373 Messnaoui and Bounahmidi (2005)

& 0–9.00 298 Cherif et al. (2000)
Speciation 0–18.0 298 Yang et al. (2016)
0–7.00 298–540 Holmes and Mesmer (1999)
0–4.37 283–353 Preston and Adams (1979)
H2 SO4 − H2 O Osmotic coefficient 0–6.00 298–328 Sippola and Taskinen (2014)
& 0–4.00 298 Messnaoui and Bounahmidi (2006)
Speciation 0–6.00 298 Christov and Moller (2004)
0–23.5 298 Irish and Chen (1971)
0–6.00 298–373 Wang et al. (2013)
0–6.00 298–473 Sippola (2012)
CaSO4 : 2H2 O − H2 SO4 Gypsum solubility 0–4.0 298–363 Wang et al. (2013)
0-0.49 298–363 Azimi and Papangelakis (2010)
0-0.79 303–353 Dutrizac (2002)
0–4.70 298–333 Shen et al. (2020)
0–6.84 283–368 Zdanovskii and Vlasov (1968)
CaSO4 − H2 SO4 Anidrhyte solubility 0–4.65 298–363 Wang et al. (2013)
0-0.79 363–378 Dutrizac (2002)
0–1.20 398-368 Shen et al. (2020)
0–4.37 283–323 Zdanovskii and Vlasov (1968)
HF − H2 O Water activity 0–6.00 298–353 EL Guendouzi et al. (2019)
0–6.00 298 Braddy et al. (1994)
H2 SiF6 − H2 O Water activity 0–3.00 298–353 El Guendouzi et al. (2015)

at different temperatures, (ii ) a pycnometer of known volume and

 β (1 )  √ √ mass is used to measure the density of the samples, (iii) the total
λi, j (I ) = i,2j −1 + (1 + α I + α I )e−α I (22)
α I concentrations of the ions and molecules present in the acid solu-
tions are calculated using a conductivity based method developed
where: by McCleskey et al. (2012) for the case of sulfuric acid, and a pH
ci,(0j) based method developed by Elmore et al. (1965) for the case of
βi,(0j ) (T ) = ai,(0j) + b(i,0j) .(T − T0 ) + (23) phosphoric acid.
T − T0
For all the experiments, 98w% sulfuric acid and 65w% phos-
phoric acid were used to prepare several samples. The tempera-
ci,(1j)
βi,(1j ) (T ) = ai,(1j) + b(i,1j) .(T − T0 ) + (24) ture and the concentration of the samples varied respectively from
T − T0 298K and 353K and from infinite dilution to 4moles/kgw for sulfu-
ric acid and from infinite dilution to 16moles/kgw for phosphoric
ci,(2j,k
)
acid.
ψi, j,k (T ) = ai,(2j,k) + b(i,2j,k) .(T − T0 ) + (25)
T − T0
Therefore, the Pitzer model involves several binary and ternary 4.2. Estimability and identification of the model parameters
interaction parameters that must be identified from experimental
measurements. However, an important question that always arises The estimability of the unknown parameters of Pitzer’s model is
is to know whether the available data allow to identify all the carried out using our recently developed method (Bouchkira et al.,
model unknown parameters or only a subset of them. To answer 2021c) based on the matrix of global sensitivities. The most es-
this question, a database of experimental measurements as well as timable parameters are then identified from the available mea-
an estimability analysis of the unknown parameters are presented surements using an optimization method based on the branch-
in the next section. and-reduce algorithm developed by Belotti et al. (2009); Sahini-
dis (2002), to minimize the following mean-square-error objective
4. Experimental database and parameter identification function to global optimality:
2
The identification of the unknown parameters of the thermo- mse
 kexp
m,T
− kmod
m,T
J = (26)
dynamic model requires experimental measurements of the most k T m
ωk
important variables. In this work, experimental measurements con-
cerning sulfuric and phosphoric acid systems have been carried out k corresponds to the measured variables, i.e., water activity, com-
in the laboratory, whereas data concerning other thermodynamic ponent molalities and solubility, ωi are weighting factors, T and
systems were collected from the literature (Table 2). m refer to the temperature and total concentrations. The accu-
racy of identification is assessed by means of confidence intervals
4.1. Laboratory measurements (Langman, 1986). Assuming that the objective function can be ex-
pressed in the following form:
Different series of experimental measurements of sulfuric and F (P ) = θ (P )T .W.θ (P ) (27)
phosphoric acid solutions have been carried out. They consist of
pH, density, electrical conductivity and molalities of the ionic and θ (P ) is the residual vector, i.e., the difference between the model
molecular species considered. These measurements were carried predictions and the measured values of the outputs, W is a weight-
out as follows: (i ) the pH and conductivity are measured by means ing matrix (equal to the identity, since we assume that all mea-
of a glass electrode multi-meter, calibrated using buffer solutions surements have the same weight in the objective function), and P

5
I. Bouchkira, A.M. Latifi, L. Khamar et al. Computers and Chemical Engineering 156 (2022) 107536

Table 3
Interaction parameters along with their confidence intervals. (A): sulfur, (B): phosphorus acid, (C): fluorine, (D): silicon, (E): calcium.
Parameters without confidence intervals are not estimable and their values are taken from the literature, or from previous works
or set equal to zero.

a C.I b C.I c C.I

(A)
(0 )
βHSO −
,SO24−
-2.471 [-3.00;-1.934] 0.007 [0.004;0.010] 0.01 -
4

βH(0SO) − ,H+ 1.046 [0.518;1.573] -0.002 [-0.003;-0.001] 0.01 -

4
(0 )
βSO 2−
,H +
-8.001 [-9.665;-6.336] 0.01 - 1684.0 [1221.1;2147.04]
4
(1 )
βHSO− ,SO2− -254.8 [-312.8;-196.8] 0.405 [0.341;0.468] 3.9.104 [3.1.104 ;4.8.104 ]
4 4

βH(1SO) − ,H+ -296.9 [-489.1;-104.8] 0.453 [0.341;0.564] 4.8.104 [2.7.104 ;7.0.104 ]

4
(1 )
βSO 2−
,H +
8.039 [4.879;11.19] -0.023 [-0.032;-0.013] 0.01 -
4
ψH+ ,HSO−4 ,SO24− 1.696 [1.142;2.249] -0.003 [-0.004;-0.001] -277.209 [-335.9;-218.4]
(B)
βH(0+),H3 PO4 4.211 [3.75;4.67] 0.134 [0.12;0.15] 2.49.103 [2.21.103 ;2.76.103 ]
βH(0+),H2 PO− -2.965 [-3.32;-2.61] -1.005 [-1.13;-0.88] 203.473 [179.06;227.89]
4

βH(0+),H5 P2 O− -186.465 [-242.4;-130.5] 2.203 [1.54;2.86] 0.198 [0.14;0.26]

8
(0 )
βH3 PO4 ,H2 PO− 0.183 [0.17;0.20] -0.088 [-0.10;-0.08] -1.25.104 [-1.37.104 ;-1.14.104 ]
4

βH(03 )PO4 ,H5 P2 O− 4.739 [4.17;5.31] -0.161 [-0.18;-0.14] 0.925 [0.81;1.04]

8

βH(02 )PO− ,H5 P2 O− 168.178 [126.1;210.2] 1.941 [1.46;2.43] 4.49.105 [3.36.105 ;5.61.105 ]
4 8
(1 )
βH+ ,H3 PO4 -2.495 [-3.24;-1.75] -0.311 [-0.40;-0.22] -4.96.104 [-6.45.104 ;-3.47.104 ]
βH(1+),H2 PO− 16.109 [12.4;19.8] 1.404 [1.08;1.73] - -
4

βH(1+),H5 P2 O− 56.786 [45.4;68.1] - - 9.14.104 [7.31.104 ;1.09.104 ]

8
(1 )
βH3 PO4 ,H2 PO− -8.458 [-9.30;-7.61] 0.064 [0.06;0.07] -1871.020 [-2058.12;-1683.92]
4

βH(13 )PO4 ,H5 P2 O− -3.712 [-4.19;-3.23] - - -3.05.104 [-3.44.104 ;-2.65.104 ]

8

βH(02 )PO− ,H5 P2 O− -4.093 [-4.91;-3.27] 0.236 [0.19;0.28] - -

4 8
ψH3 PO4 ,H3 PO4 ,H3 PO4 6.826 [5.46;8.19] - - - -
βH(03 )PO4 ,H3 PO4 - - - - -134.046 [-135.39;-132.71]
(C)
βH(0+),HF - - -0.007 [-0.0078;-0.0062] - -
βH(0F)− ,HF 0.208 [0.19;0.23] - - - -
2

βF(0−),HF - - 0.013 [0.0118;0.0142] - -

βH(0F)− ,HF -0.392 [-0.44;-0.35] . - - -
2
ψHF,HF,HF -3.081 [-4.01;-2.16] - - - -
βH(0F,H
)
F
- - - - -84.628 [-93.94;-75.32]
(D)
βH(0+),H2 SiF6 - - 0.072 [0.064;0.079] - -
βH(0+),SiF 2− - - -1.338 [-1.45;-1.21] -9.071.104 [-1.00.105 ;-8.07.104 ]
6

βH(0)SiF ,SiF 2− 12.112 [10.7;13.4] 0.077 [0.067;0.086] -3.746 [-4.08;-3.40]

2 6 6

βH(02 )SiF6 ,H2 SiF6 - - 0.005 [0.004;0.0055] - -

βH(1+),H2 SiF6 -4.602 [-5.01;-4.18] -7.095 [-8.15;-6.030] -7.728.105 [-8.65.105 ;-6.80.105 ]
βH(1+),H2 SiF6 4.901 [4.31;5.48] - - -9.676.104 [-1.11.104 ;-8.22.104 ]
βH(1+),SiF 2− -63.290 [-72.7;-53.7] - - -8.341.104 [-9.42.104 ;-7.25.104 ]
6

βH(1)SiF ,SiF 2− -0.149 [-0.16;-0.12] - - - -

2 6 6
(1 )
βH2 SiF6 ,H2 SiF6 0.107.10−6 [0.01.10−6 ;2.78.10−6 ] 0.040 [0.036;0.0436] - -
(E)
(0 )
βHSO −
,CaSO4
0.920 [0.81;1.02] -0.030 [-0.033;-0.026] -3272.216 [-3632.1;-2912.2]
4
(0 )
βHSO −
,Ca2+
1.010 [0.91;1.10] 0.026 [0.023;0.028] - -
4
(0 )
βSO 2−
,CaSO4
-1.801 [-2.01;-1.58] 0.010 [0.0088;0.0112] - -
4
(0 )
βSO 2−
,Ca2+
-1.037 [-1.19;-0.88] 0.066 [0.0561;0.0759] -1075.927 [-1194.2;-957.5]
4

βH(0+),CaSO4 0.697 [0.60;0.78] 0.033 [0.028;0.037] - -

(1 )
βSO 2−
,Ca2+
- - - - 23487.762 [20904.1;26071.4]
4

βH(1+),CaSO4 -7.917 [-8.78;-7.04] - - - -

βH(1+),Ca2+ - - - - -912.578 [-1022.0;-803.0]
βC(aSO
1)
4 ,C a
2+ -26.026 [-28.8;-23.1] - - - -
ψCaSO4 ,CaSO4 ,CaSO4 0.689 [0.61;0.76] - - - -

is the vector of the unknown parameters. Assuming that the mea- parameters, and J is the Jacobian matrix of the vector θ (P ). This
surement errors are independent and normally distributed, the co- approximation is more accurate when the non-linearities are not
variance matrix C of the least square problem is approximated as strong.
follows (Langman, 1986): The Jacobian matrix is computed using a local ”One-At-Time”
F (P ∗ ) T −1 (OAT) method. It consists in perturbing the value of each pa-
C≈ (J J ) (28) rameter Pj∗ by 10% forward and backward, then, the centered fi-
d−n
nite difference method is used to approximate the elements of
P ∗ is the vector of parameters that minimizes F (P ), n and d are re-
J. It is worth noting that 10% disturbance is chosen based on
spectively the number of experiments and the number of unknown

6
I. Bouchkira, A.M. Latifi, L. Khamar et al. Computers and Chemical Engineering 156 (2022) 107536

many previous works that have addressed the same problems (Wu the least estimable ones. To answer this question, we have devel-
et al., 2011; Shahmohammadi and McAuley, 2018; Yao et al., 2003; oped, in one of our recent works (Bouchkira et al., 2021c), a new
Fysikopoulos et al., 2019; Benyahia et al., 2013; Ngo et al., 2014; global sensitivity-based estimability analysis approach. The latter
Lei et al., 2013; Cardenas et al., 2021, 2020). On the other hand, J was then used to determine the number of optimal parameters to
is used to analyze whether performing additional experiments will be identified for each considered thermodynamic subsystem (i.e.,
increase the number of parameters they can estimate (estimabil- sulfur, phosphorus, fluorine, silica and calcium), based on the cal-
ity analysis), while the Fisher Information Matrix (F IM ) (expressed culation of an estimability threshold.
as F IM = J T V −1 J, where V is the variance-covariance matrix of the In the case of interaction parameters based on sulfur, the cal-
outputs) is used to find out how these proposed experiments will culated threshold of estimability is equal to 0.05, which means
improve the precision of the estimates of the important parame- that for a parameter to be estimable, a variation of 10% of its
ters of the model (optimal design of experiments). Furthermore, value produces a variation of at least 2.24% in model outputs. With
F IM is often utilized to determine the parameter confidence inter- this threshold value, 17 interaction parameters can be estimated
vals/regions (see Section 4.2). from the available measurements. In the other cases, the calculated
⎛ ∂θ1 ∂θ1 ⎞ threshold of estimability, the minimum variation of the model out-
∂ P1 ··· ∂ Pd puts produced by a perturbation of 10% of a parameter, and the
⎜ ⎟
J = ⎝ ... ..
.
..
. ⎠ number of estimable parameters are respectively 0.02, 1.41%, and
∂θn ∂θn 34 for the phosphorus; 0.013, 1.14% and 6 for fluorine; 0.015, 1.34%
∂ P1 ··· ∂ Pd and 17 for silica; and 0.053, 1.34% and 17 for calcium.
The uncertainty on a parameter j is determined by means of As expected, the estimability analysis results show that the
Eq. (29), where c j j is the jth diagonal element of the covariance available data do not allow to identify all the unknown parame-
matrix C. The latter is generally estimated from experimental mea- ters of the developed thermodynamic model. Indeed, the unknown
surements and then used to estimate confidence intervals. How- parameters of the Pitzer model describe binary and ternary inter-
ever, since in this work the estimation of C was not easy to achieve, actions between different ionic or molecular species in the aque-
it is approximated using F (P ∗ ) which represents the residues be- ous phase. The non-estimability of a parameter involved in the de-
tween the model predictions and the experiments used. Moreover, scription of an interaction between species, means that the used
the division by P ∗ allows us to quantify the precision of the iden- experimental measurements do not contain enough information
tification of a parameter P . about that interaction. One possible improvement is to design op-
√ timal experiments that are more sensitive to interactions whose
c j j t1−α /2,ν
Pj = ± .100% (29) model involves the least estimable parameters. As an example,
Pj∗ pH measurements are usually used to improve the estimability of
t1−α /2,ν is deduced from the Student table, it corresponds to unknown parameters involving interactions between H + ions and
the probability of 1 − α /2 that the true value of the parameter is other species. Finally, the values of the unknown parameters which
within the confidence interval. ν is the degree of freedom given by cannot be estimated from the available measurements are either
the difference between the number of data points and the number fixed from previous works or from literature, or simply set equal
of unknown parameters, i.e., v = n − d. The confidence intervals are to zero.
given by: The parameter identification problems are then implemented
    within GAMS environment and solved using both a Dell Precision
Pj ∈ Pj∗ − c j j t1−α /2,ν ; Pj∗ + c j j t1−α /2,ν (30) T7810 Bi-Xeon 12 x Core 64GB work station, and the NEOS Server
On the other hand, the Pearson-product-moment-coefficient is cal- which provides access to more than sixty state-of-the-art solvers
culated to evaluate the accuracy and reliability of the predictions in more than a dozen optimization categories, hosted by machines
of the developed thermodynamic-based model as (Keith, 2014): at Arizona State University, the University of Klagenfurt in Austria,
    and the University of Minho in Portugal (Czyzyk et al., 1999; Dolan,
i, j Simj − Simj kei j − kei j 2001; Gropp and Moré, 1997). On the other hand, Baron optimizer,
r=  2  2 (31) based on the branch-and-reduce algorithm, is used to solve the

i, j (Si j − Si j ) i, j (ki j − ki j )
m m e e
problems to global optimality. The CPU time for each minimization
problem is set to about three hours. The number of parameters in-
where Simj is the model prediction, kei j is the corresponding experi-
volved in each minimization problem can be found in Tables 3A-E.
mental value, Simj and kei j are respectively their average over the to- The values of the identified parameters are listed in Table 3 along
tal number of available measurements performed at different mo- with their confidence intervals. It is worth mentioning that param-
lalities and temperatures. The indices i and j refer to molalities and eters without interval confidence are not estimable, and their val-
temperatures respectively. ues are taken from previous works or from the literature. These
parameters are then used to compare the predictions of the devel-
5. Results and discussions oped model to the available measurements.
The results are presented in Figs. 2A-H, 3 A-H and 4 A-H. To
5.1. Estimability analysis and parameter identification assess the accuracy of the model predictions, the Pearson product-
moment coefficient r is also computed in each case. Its high values
In many previous studies which deal with predictive models highlight the accuracy and the reliability of the model.
that involve unknown parameters to be identified from available The developed model is then exploited in the multi-objective
experiments, the authors assume that the available data contain optimization problem to predict the values of the objective func-
the necessary information to accurately identify them. However, tions. The  -constraint method (Deb, 2014) is used to transform
it is well recognized that it is not always the case. Indeed, due the multi-objective optimization problem to a single-objective op-
to limited data, noisy measurements, or sometimes correlated de- timization problem which is then solved several times accurately
signs of experiments, the estimation of the unknown parame- by means of a global optimization solver to determine the Pareto
ters may not be accurate and the underlying question is to know front of optimal solutions. The multi-objective optimization prob-
which parameters are estimable from the available measurements, lem is also implemented within GAMS environment and solved by
and eventually to design appropriate experiments to determine means of Baron optimizer on the same computers described above.

7
I. Bouchkira, A.M. Latifi, L. Khamar et al. Computers and Chemical Engineering 156 (2022) 107536

Fig. 2. Model vs measurements. (A-D): speciation of sulfuric acid. (E-H): speciation of phosphoric acid.

8
I. Bouchkira, A.M. Latifi, L. Khamar et al. Computers and Chemical Engineering 156 (2022) 107536

Fig. 3. Model vs measurements. (A-D):gypsum solubility in sulfuric acid. (E-H): anhydrite solubility in sulfuric acid.

9
I. Bouchkira, A.M. Latifi, L. Khamar et al. Computers and Chemical Engineering 156 (2022) 107536

Fig. 4. Model vs measurements. (A-H): water in different systems. (1: H2 SO4 − H2 O, 2: H3 PO4 − H2 O, 3: H2 SiF6 − H2 O, 4: HF − H2 O).

10
I. Bouchkira, A.M. Latifi, L. Khamar et al. Computers and Chemical Engineering 156 (2022) 107536

Table 4
Four solutions from the Pareto front.

w1 w2 w3 T opt (◦C ) Sc ( % ) Qcool (kJ/hr )

1 0.73 0.27 0 76.3 2.12 -1.333. 105

2 0.70 0.30 0 77.1 2.18 -1.351. 105
3 0.67 0.33 0 77.6 2.23 -1.355. 105
4 0.63 0.37 0 78.2 2.54 -1.358. 105

Furthermore, the excess of sulfuric acid (with respect to the

stoichiometric amount) lies between 2% and 3%, which is in good
agreement on the one hand with the current operating condi-
tions of the industrial units, and on the other hand with the con-
ditions favorable to the crystallization of gypsum described by
Becker (1989). It is noteworthy that the optimal solutions for sul-
furic acid distributions are not uniform. In fact, as the conversion
Fig. 5. Pareto front of the optimization problem with four points enumerated 1–4
and reported in Table 4. rate of the digestion increases from the first to the last section, the
demand for sulfuric acid varies in the opposite direction and de-
creases from the first to the third section.
The CPU time is set to about seventeen hours since the number of Optimal temperatures are also reported in Table 4, their values
optimal solutions which constitute the Pareto front is equal to 50, allow to determine the overall heat quantity to be removed via the
meaning that the multi-objective optimization problem is trans- flash cooler to ensure good conditions for gypsum crystallization
formed into 50 single optimization problems with 20 min CPU and also to limit the formation of bassanite which undermines the
time each. quality of the filtration unit downstream of the digestion tank.

5.3. Decision-making aid method

5.2. Optimization results
Once the Pareto front has been determined, the goal is to rank
The multi-objective optimization problem is formulated as the
all the optimal solutions in order to select the best one to im-
following single-objective optimization problem:
plement. A multi-objective decision-making aid method is there-
fore needed. In this work, the multi attribute-utility theory is used
Maxx − f 1 (x )
(Mateo, 2012). It is based on the calculation of individual utility
subject to f 2 (x ) ≥  functions u1 and u2 , which mainly scale the objective functions f1
g j (x ) = 0 j ∈ [1; J] and f2 to the same magnitude, and then to rank the optimal solu-
hk ( x ) ≥ 0 k ∈ [1; K] tions through the calculation of the multi-attribute utility function
U.
xi(L ) ≤ xi ≤ xi(U ) i ∈ [1; N] α1
f1max − f1 (x )
where f1 and f2 are the two conflicting objectives predicted us- u1 ( x ) = (32)
f1max − f1min
ing the developed thermodynamic model, x is the vector of deci-
α2
sion variables, i.e., [w1 , w2 , w3 , Sc , Qcool ], x(L ) and x(U ) their lower f2 (x ) − f2min
and upper bounds. g = (g1 , g2 , . . . , gJ )T and J are the equality con- u2 ( x ) = (33)
f2max − f2min
straints and their number respectively, h = (h1 , h2 , . . . , hK )T and K
refer to the inequality constraints and their number. f1max and f2max are the maximum values of the objectives of the
The epsilon-constraint optimization problem is then imple- Pareto front and f1min and f2min are their minimum values. α1 and
mented and solved within GAMS environment. Moreover, a branch- α2 are related to the relative tolerance on the individual utility
and-reduce algorithm (Belotti et al., 2009; Sahinidis, 2002) is used functions. The multi-attribute utility function is defined as:
to solve the problem to global optimality. It should be noted that
U (x ) = w1 .u1 (x ) + w2 .u2 (x ) (34)
the method has the advantage of finding any Pareto optimal solu-
tion, even for non-convex problems. However, it is not always easy where w1 and w2 are usually set by a decision-maker according
to set the lower limit  so that the problem has feasible solutions. to his/her preferences. They express the weight of each objective
The Pareto front computed is presented in Fig. 5 and shows that in the decision. The values of the multi-attribute utility function
(i) the saturation index of brushite (syncristallized losses) is nega- allow then to rank the optimal solutions and to select the best one
tive for all the optimal solutions meaning that the thermodynamic to be implemented.
conditions are not favorable for its production, (ii) the saturation Figure 6A-B show the effect of the relative tolerances α1 and
index of gypsum is positive with higher values for all the optimal α2 on the individual utility functions. They represent the sensitiv-
solutions meaning that the thermodynamic conditions are very fa- ity of u1 and u2 with respect to f1 and f2 . Their values are set
vorable for its production and consequently for increasing the pro- in our case to α1 = 0.5 and α2 = 0.5. Meanwhile, the values of the
ductivity of the industrial unit. weighting parameters w1 and w2 are taken equal to 0.5 to consider
On the other hand, solutions are taken from the Pareto front that the two objective functions have the same importance in the
and detailed in Table 4. They show that the optimal distribution optimization problem.
of sulfuric acid in the three sections varies between 63% and 73% The individual utility function U (x ) was calculated for each
(w1 ) for the first section, between 27% and 37% (w2 ) for the sec- point of the Pareto front and used to rank the optimal solutions.
ond and 0% (w3 ) for the last section. These values are very con- The first two best solutions are shown in Table 5 and in Fig. 7,
sistent with the current distributions used on industrial scale, in along with industrial practices from (Bouchkira and El Fariq, 2017;
particular within a phosphoric acid plant located in Jorf Lasfar in El Bouzidi, 2016) and other optimal solutions from the Pareto front.
Morocco (Bouchkira and El Fariq, 2017; El Bouzidi, 2016). It is important to point out that the total losses seem very small

11
I. Bouchkira, A.M. Latifi, L. Khamar et al. Computers and Chemical Engineering 156 (2022) 107536

Fig. 6. Relative tolerance factors αi effect. (a) Saturation index of gypsum, (b) Saturation index of Brushite.

Fig. 7. Comparison between optimal solution, industrial practices and two best solutions.

Table 5 5.4. Surrogate model

Industrial measurements and optimal solutions.

w1 w2 w3 SL(% ) UL(% ) Reference A surrogate model is developed in order to compare the opti-
1 0.50 0.20 0.30 0.94 0.06 El Bouzidi (2016) mal solutions with the current practice at the plant. Indeed, this
2 0.40 0.20 0.40 0.73 0.12 comparison cannot be made directly through the values of the ob-
3 0.40 0.30 0.30 0.87 0.12 jective functions used in optimization, since it is practically impos-
4 0.70 0.30 0 0.62 0.14 Bouchkira and El Fariq (2017) sible to calculate the saturation indices from the industrial prac-
5 0.64 0.36 0 0.70 0.15
tice.
6 0.75 0.25 0 0.64 0.12
7 0.75 0.17 0.08 0.68 0.14 For this reason, a fractional design of experiments (Park, 2007)
8 0.60 0.30 0.10 0.63 0.12 is used to develop the surrogate model, which enables to pre-
9 0.65 0.35 0 0.64 0.12 This work dict the syncrystallized as well as the unattacked losses of phos-
10 0.63 0.37 0 0.62 0.14
phate. Given the complexity of the process and lack of indus-
trial measurements, two factors (i.e., temperature and excess sul-
furic acid) are assumed to be known. Therefore, only three fac-
tors out of five, i.e., w1 , w2 , w3 , are considered with two
and one might be tempted to ignore them. However, the industrial
levels.
phosphoric acid plant located in Jorf Lasfar in Morocco (studied in
The design of experiments factors (i.e., w1 , w2 and w3 ) along
this work), includes 16 units treating 1400 tons of phosphate ore
with their levels (i.e., L(−1 ) and L(+1 ) are reported in Table 6. The
per day each. Therefore, 1% loss corresponds to 224 tons of phos-
phate loss per day, which are significant and cannot be neglected.

12
I. Bouchkira, A.M. Latifi, L. Khamar et al. Computers and Chemical Engineering 156 (2022) 107536

Table 6 Declaration of Competing Interest

Design of experiments levels: L(+1) and L(−1) . SL and UL refer to
syncristallized and unattacked losses.
The authors declare that they have no known competing finan-
w1 w2 w3 SL(% ) UL(% ) cial interests or personal relationships that could have appeared to
1 +1 +1 -1 0.62 0.14 influence the work reported in this paper.
2 +1 -1 -1 0.64 0.12
3 -1 -1 +1 0.94 0.06 CRediT authorship contribution statement
4 -1 +1 +1 0.87 0.12
L(+1) [0.5;0.7] [0.2;0.3] [0.2;0.4] - -
L(−1) [0.4;0.5] [0.1;0.2] [0; 0.2] - - Ilias Bouchkira: Conceptualization, Data curation, Soft-
ware, Writing – original draft, Writing – review & editing.
Abderrazak M. Latifi: Conceptualization, Funding acquisition,
Methodology, Supervision, Writing – review & editing. Lhachmi
surrogate model adopted is expressed as :
Khamar: Investigation, Writing – review & editing. Saad Benjel-
  loun: Resources, Writing – review & editing.
Y mod = a0 + bi xi + ci, j xi x j (35)
i i j
Acknowledgments

θ = (X T . X )−1 . X T .Y exp (36) The authors are grateful to OCP Group for its financial support.
where Y mod refers to the surrogate model predictions. Y exp is the
Supplementary material
vector of measured values. X is the matrix of the experimental de-
sign. xi and x j are the levels of the three considered variables. θ
Supplementary material associated with this article can be
is the vector of the unknown coefficients of the surrogate model
found, in the online version, at doi10.1016/[Link].2021.
given as:
107536.
θ = [a0 b1 b2 b3 c1,2 c1,3 c2,3 ]t (37)
References
Figure 7 shows a comparison between industrial losses taken
Abu-Eishah, S.I., Abu-Jabal, N.M., 2001. Parametric study on the production of phos-
from (Bouchkira and El Fariq, 2017; El Bouzidi, 2016), the opti- phoric acid by the dihydrate process. Chem. Eng. J. 81 (1–3), 231–250. doi:10.
mal solutions from the Pareto front and the first two best opti- 1016/S1385-8947(0 0)0 0166-2.
mal solutions. The computed solutions are very close to the in- Al-Fariss, T., Razik, S.A., Ghasem, N.M., Abdelaleem, F., Elnashaie, S.S., 1993. Com-
puter program for the mass and heat balance calculations for the wet process
dustrial measurements which are obtained by performing sev-
phosphoric acid (hemi-hydrate) and its applications to Saudi phosphate rock.
eral run-to-run experiments on the plant by changing the val- Fert. Res. 35 (3), 169–175. doi:10.10 07/BF0 0750635.
ues of the decision variables from one run to the next to im- Al-Thyabat, S., Zhang, P., 2015. In-line extraction of REE from Dihydrate (DH) and
HemiDihydrate (HDH) wet processes. Hydrometallurgy 153, 30–37. doi:10.1016/
prove the productivity and minimize the phosphate losses. More-
[Link].2015.01.010.
over, the values of the decision variables used to carry out in- Alfariss, T.F., Özbelge, H., El-Shall, H., 1992. Process technology for phosphoric acid
dustrial experiments, are sub-optimal and must be adjusted de- production in Saudi Arabia. J. King Saud Univ.-[Link]. 4 (2), 239–254. doi:10.
pending on the composition of the phosphate ore which further 1016/S1018-3639(18)30567-1.
Amin, M., Ali, M., Kamal, H., Youssef, A., Akl, M., 2010. Recovery of high grade phos-
increases the number of runs. It is easy to understand under phoric acid from wet process acid by solvent extraction with aliphatic alcohols.
these conditions that obtaining these measurements is very costly Hydrometallurgy 105 (1–2), 115–119. doi:10.1016/[Link].2010.08.007.
since it involves significant labor cost, consumes time and huge Azimi, G., Papangelakis, V., 2010. The solubility of gypsum and anhydrite in simu-
lated laterite pressure acid leach solutions up to 250 ◦ C. Hydrometallurgy 102
quantities of reagents, decreases the productivity of the plant and (1–4), 1–13. doi:10.1016/[Link].20 09.12.0 09.
poses some safety issues. Modeling and multi-criteria optimization Ball, J.W., Nordstrom, D.K., 1991. WATEQ4F–User’s Manual with Revised Thermody-
have proven to be effective since they provide optimal solutions namic Data Base and Test Cases for Calculating Speciation of Major, Trace and
Redox Elements in Natural Waters. Technical Report.
that minimize the number, duration, cost of runs, and security Becker, P., 1989. Phosphates and Phosphoric Acid: Raw Materials, Technology, and
risks. Economics of the Wet Process. Revised and Expanded, vol. 6. Marcel Dekker,
Inc..
Belboom, S., Szöcs, C., Léonard, A., 2015. Environmental impacts of phosphoric acid
6. Conclusion production using di-hemihydrate process: a Belgian case study. J. Clean. Prod.
108, 978–986. doi:10.1016/[Link].2015.06.141.
Belotti, P., Lee, J., Liberti, L., Margot, F., Wächter, A., 2009. Branching and bounds
Optimal operating conditions of the digestion tank of an in- tighteningtechniques for non-convex MINLP. Optim. Methods Softw. 24 (4–5),
dustrial phosphoric acid production process are determined by 597–634. doi:10.1080/10556780903087124.
means of a multi-objective optimization method. The objectives Benyahia, B., Latifi, M.A., Fonteix, C., Pla, F., 2013. Emulsion copolymerization of
styrene and butyl acrylate in the presence of a chain transfer agent. Part 2:
used (i.e., saturation indices) have proven to be effective in min- parameters estimability and confidence regions. Chem. Eng. Sci. 90, 110–118.
imizing the phosphate losses and increasing the acid productiv- doi:10.1016/[Link].2012.12.013.
ity. Furthermore, the two inequality constraints on the temperature Bharathi, G., Prabhakaran, D., Kannadasan, T., 2013. Analysis and simulation of di-
hydrate process for the production of phosphoric acid (reactor section). Am. J.
and on the excess of sulfuric acid above the stoichiometric amount Eng. Res 2, 1–8. doi:10.1007/BF02697694.
required by the reactions were relevant to ensure an optimal oper- Bichri, A., Kamzon, M.A., Abderafi, S., 2020. Artificial neural network to predict the
ation of the digestion tank. performance of the phosphoric acid production. Procedia Comput. Sci. 177, 444–
449. doi:10.1016/[Link].2020.10.060.
The optimization results are consistent with the measurements
Bouchkira, I., Benjelloun, S., Khamar, L., Latifi, A.M., 2021. Thermodynamic-based
carried out on an industrial plant, thus showing that optimiza- model for the prediction of the fouling phenomena in a wet phosphoric acid
tion is a powerful tool for reducing the cost and duration of process. Chem. Eng. Trans. 86, 1273–1278. doi:10.3303/CET2186213.
Bouchkira, I., El Fariq, Y., 2017. Improvement of the 28% phosphoric acid production
industrial experiments under maximum safety conditions. How-
units. MSc thesis ENSA, Khouribga, Morocco.
ever, it is noteworthy that some side reactions can occur due to Bouchkira, I., Latifi, A., Khamar, L., Benjelloun, S., 2021. Process modeling and multi-
the impurities present in the phosphate ore which are not taken criteria optimization of an industrial phosphoric acid wet-process. Comput.
into account in this work. Their consideration in future works Aided Chem. Eng. 50, 499–504. doi:10.1016/B978- 0- 323- 88506- 5.50079- 6.
Bouchkira, I., Latifi, A.M., Khamar, L., Benjelloun, S., 2021. Global sensitivity based
would undoubtedly improve the performance of the digestion estimability analysis for the parameter identification of Pitzer’s thermodynamic
tank. model. Reliab. Eng. Syst. Saf. 207, 107263. doi:10.1016/[Link].2020.107263.

13
I. Bouchkira, A.M. Latifi, L. Khamar et al. Computers and Chemical Engineering 156 (2022) 107536

Bouchkira, I., Latifi, A.M., Khamar, L., Benjelloun, S., 2021. Multi-objective optimiza- Keith, T.Z., 2014. Multiple Regression and Beyond: An Introduction to Multiple Re-
tion of the digestion tank of an industrial phosphoric acid manufacturing pro- gression and Structural Equation Modeling. Routledge.
cess. In: 2020 6th IEEE Congress on Information Science and Technology (CiSt). Lampila, L.E., 2013. Applications and functions of food-grade phosphates. Ann. N. Y.
IEEE, pp. 389–394. doi:10.1109/CiSt49399.2021.9357206. Acad. Sci. 1301 (1), 37–44.
Braddy, R., McTigue, P., Verity, B., 1994. Equilibria in moderately concentrated aque- Langman, M., 1986. Towards estimation and confidence intervals. Br. Med. J. (Clin.
ous hydrogen fluoride solutions. J. Fluor. Chem. 66 (1), 63–67. doi:10.1016/ Res. Ed) 292 (6522), 716. doi:10.1136/bmj.292.6522.716.
0022- 1139(93)02895- L. Lei, M., Vallières, C., Grévillot, G., Latifi, M.A., 2013. Thermal swing adsorption pro-
Cardenas, C., Latifi, A. M., Vallières, C., Marsteau, S., Sigot, L., 2021. Analysis of an cess for carbon dioxide capture and recovery: modeling, simulation, parame-
industrial adsorption process based on ammonia chemisorption: modeling and ters estimability, and identification. Ind. Eng. Chem. Res. 52 (22), 7526–7533.
simulation. 10.1016/[Link].2021.107474. doi:10.1021/ie3029152.
Cardenas, C., Marsteau, S., Sigot, L., Vallières, C., Latifi, A.M., 2020. Analysis of an Leikam, D.F., Achorn, F.P., 2005. Phosphate fertilizers: production, characteris-
industrial adsorption process based on ammonia chemisorption: modeling and tics, and technologies. Phosphorus Agric. Environ. 46, 23–50. doi:10.2134/
simulation. In: Computer Aided Chemical Engineering, vol. 48. Elsevier, pp. 625– agronmonogr46.c2.
630. doi:10.1016/B978- 0- 12- 823377- 1.50105- 1. Li, R., Wang, J.J., Zhou, B., Awasthi, M.K., Ali, A., Zhang, Z., Lahori, A.H., Ma-
Cherif, M., Mgaidi, A., Ammar, M.N., Abderrabba, M., Fürst, W., 20 0 0. Modelling of har, A., 2016. Recovery of phosphate from aqueous solution by magnesium
the equilibrium properties of the system H3PO4–H2O: representation of VLE oxide decorated magnetic biochar and its potential as phosphate-based fertilizer
and liquid phase composition. Fluid Phase Equilib. 175 (1–2), 197–212. doi:10. substitute. Bioresour. Technol. 215, 209–214. doi:10.1016/[Link].2016.02.125.
1016/S0378-3812(0 0)0 0458-1. Ma, H., Feng, X., Deng, C., 2018. Water–phosphorus nexus for wet-process phospho-
Chidambaram, S., Karmegam, U., Sasidhar, P., Prasanna, M.V., Manivannan, R., ric acid production. Ind. Eng. Chem. Res. 57 (20), 6968–6979. doi:10.1021/acs.
Arunachalam, S., Manikandan, S., Anandhan, P., 2011. Significance of saturation iecr.7b05399.
index of certain clay minerals in shallow coastal groundwater, in and around Mateo, J.R.S.C., 2012. Multi-attribute utility theory. In: Multi Criteria Analysis in the
Kalpakkam, Tamil Nadu, India. J. Earth Syst. Sci. 120 (5), 897–909. Renewable Energy Industry. Springer, pp. 63–72.
Cho, H.J., Yeo, Y.-K., Park, W.H., Moon, B.K., 1996. Modeling and simulation of a wet McCleskey, R.B., Nordstrom, D.K., Ryan, J.N., Ball, J.W., 2012. A new method of cal-
hemihydrate phosphoric acid process. Korean J. Chem. Eng. 13 (6), 585–595. culating electrical conductivity with applications to natural waters. Geochim.
doi:10.1007/BF02706025. Cosmochim. Acta 77, 369–382. doi:10.1016/[Link].2011.10.031.
Christov, C., Moller, N., 2004. Chemical equilibrium model of solution behavior and Messnaoui, B., Bounahmidi, T., 2005. Modeling of excess properties and vapor–liquid
solubility in the H-Na-K-OH-Cl-HsO4-SO4-H2O system to high concentration equilibrium of the system H3PO4–H2O. Fluid Phase Equilib. 237 (1–2), 77–85.
and temperature. Geochim. Cosmochim. Acta 68 (6), 1309–1331. doi:10.1016/j. doi:10.1016/j.fluid.20 05.08.0 02.
gca.2003.08.017. Messnaoui, B., Bounahmidi, T., 2006. On the modeling of calcium sulfate solubility
Czyzyk, J., Wisniewski, T., Wright, S.J., 1999. Optimization case studies in the NEOS in aqueous solutions. Fluid Phase Equilib. 244 (2), 117–127. doi:10.1016/j.fluid.
guide. SIAM Rev. 41 (1), 148–163. doi:10.1137/S0036144598334874. 2006.03.022.
Deb, K., 2014. Multi-objective optimization. Search Methodologies. New York: Ngo, V.V., Michel, J., Gujisaite, V., Latifi, A., Simonnot, M.-O., 2014. Parameters de-
Springer. scribing nonequilibrium transport of polycyclic aromatic hydrocarbons through
Dolan, E. D., 2001. NEOS Server 4.0 Administrative Guide. arXiv preprint cs/0107034. contaminated soil columns: estimability analysis, correlation, and optimization.
Dorozhkin, S.V., 1996. Fundamentals of the wet-process phosphoric acid production. J. Contam. Hydrol. 158, 93–109. doi:10.1016/[Link].2014.01.005.
1. Kinetics and mechanism of the phosphate rock dissolution. Ind. Eng. Chem. Papadopoulos, A.I., Seferlis, P., 2009. Generic modelling, design and optimization
Res. 35 (11), 4328–4335. doi:10.1021/ie960092u. of industrial phosphoric acid production processes. Chem. Eng. Process. 48 (1),
Dutrizac, J., 2002. Calcium sulphate solubilities in simulated zinc processing 493–506. doi:10.1016/[Link].2008.06.011.
solutions. Hydrometallurgy 65 (2–3), 109–135. doi:10.1016/S0304-386X(02) Park, G.-J., 2007. Analytic Methods for Design Practice. Springer Science & Business
0 0 082-8. Media.
El Bouzidi, S., 2016. Improvement of two phosphoric acid units. MSc thesis, ENSA, Peng, Y., Zhu, Z., Braatz, R.D., Myerson, A.S., 2015. Gypsum crystallization dur-
Khouribga, Morocco. ing phosphoric acid production: modeling and experiments using the mixed-
EL Guendouzi, M., Faridi, J., Khamar, L., 2019. Chemical speciation of aqueous hy- solvent-electrolyte thermodynamic model. Ind. Eng. Chem. Res. 54 (32), 7914–
drogen fluoride at various temperatures from 298.15 K to 353.15 K. Fluid Phase. 7924. doi:10.1021/[Link].5b01763.
Equilib. 499, 112244. doi:10.1016/j.fluid.2019.112244. Pitzer, K.S., 2018. Activity Coefficients in Electrolyte Solutions. CRC press.
El Guendouzi, M., Rifai, A., Skafi, M., 2015. Properties of fluoride in wet phosphoric Preston, C.M., Adams, W., 1979. A laser raman spectroscopic study of aqueous or-
acid processes: fluorosilicic acid in an aqueous solution of H2SiF6–H20 at tem- thophosphate salts. J. Phys. Chem. 83 (7), 814–821. doi:10.1021/j100470a011.
peratures ranging from 298.15 k to 353.15 k. Fluid Phase Equilib. 396, 43–49. Rabadjieva, D., Sezanova, K., Gergulova, R., Titorenkova, R., Tepavitcharova, S., 2020.
doi:10.1016/j.fluid.2015.03.014. Precipitation and phase transformation of dicalcium phosphate dihydrate in
Elmore, K., Hatfield, J., Dunn, R., Jones, A., 1965. Dissociation of phosphoric acid electrolyte solutions of simulated body fluids: thermodynamic modeling and ki-
solutions at 25. J. Phys. Chem. 69 (10), 3520–3525. doi:10.1021/j100894a045. netic studies. J. Biomed. Mater. Res. Part A 108 (8), 1607–1616. doi:10.1002/jbm.
Fysikopoulos, D., Benyahia, B., Borsos, A., Nagy, Z.K., Rielly, C.D., 2019. A framework a.36929.
for model reliability and estimability analysis of crystallization processes with Rickard, D.A., 20 0 0. Review of phosphorus acid and its salts as fertilizer materials.
multi-impurity multi-dimensional population balance models. Comput. Chem. J. Plant Nutr. 23 (2), 161–180. doi:10.1080/01904160 0 093820 06.
Eng. 122, 275–292. doi:10.1016/[Link].2018.09.007. Sahinidis, N.V., 2002. Global optimization and constraint satisfaction: the branch-
Geeson, M.B., Cummins, C.C., 2018. Phosphoric acid as a precursor to chemicals tra- and-reduce approach. In: International Workshop on Global Optimization and
ditionally synthesized from white phosphorus. Science 359 (6382), 1383–1385. Constraint Satisfaction. Springer, pp. 1–16. doi:10.1007/978- 3- 540- 39901- 8_1.
doi:10.1126/science.aar6620. Shahmohammadi, A., McAuley, K.B., 2018. Sequential model-based a-optimal design
Gilmour, R., 2013. Phosphoric Acid: Purification, Uses, Technology, and Economics. of experiments when the fisher information matrix is noninvertible. Ind. Eng.
CRC Press. Chem. Res. 58 (3), 1244–1261. doi:10.1021/[Link].8b03047.
Gioia, F., Mura, G., Viola, A., 1977. Analysis, simulation, and optimization of the Shen, L., Sippola, H., Li, X., Lindberg, D., Taskinen, P., 2020. Thermodynamic mod-
hemihydrate process for the production of phosphoric acid from calcareous eling of calcium sulfate hydrates in a CaSO4–H2SO4–H2O system from 273.15
phosphorites. Ind. Eng. Chem. Process Des. Dev. 16 (3), 390–399. doi:10.1021/ to 473.15 K up to 5 m sulfuric acid. J. Chem. Eng. Data 65 (5), 2310–2324.
i260063a025. doi:10.1021/[Link].9b00829.
Grema, A., Imam, Y., Mohammed, H., 2018. Modeling and simulation of hemihy- Simoes, M.C., Hughes, K.J., Ingham, D.B., Ma, L., Pourkashanian, M., 2017. Tempera-
drate phosphoric acid plant. Arid Zone J. Eng. [Link]. 14 (2), 169– ture dependence of the parameters in the Pitzer equations. J. Chem. Eng. Data
171. 62 (7), 20 0 0–2013. doi:10.1021/[Link].7b0 0 022.
Gropp, W., Moré, J.J., 1997. Optimization Environments and the NEOS Server. Tech- Sippola, H., 2012. Thermodynamic modelling of concentrated sulfuric acid solutions.
nical Report. Argonne National Lab., IL (United States). Calphad 38, 168–176. doi:10.1016/[Link].2012.06.008.
Gurgul, S.J., Seng, G., Williams, G.R., 2019. A kinetic and mechanistic study into the Sippola, H., Taskinen, P., 2014. Thermodynamic properties of aqueous sulfuric acid.
transformation of calcium sulfate hemihydrate to dihydrate. J. Synchrotron Ra- J. Chem. Eng. Data 59 (8), 2389–2407. doi:10.1021/je4011147.
diat. 26 (3), 774–784. doi:10.1107/S160 05775190 01929. Slack, A., James, G., 1973. Ammonia. Fertilizer science and technology series. Marcel
Gustafsson, J.P., 2011. Visual MINTEQ 3.0 User Guide. KTH, Department of Land and Deckker Inc., New York.
Water Recources, Stockholm, Sweden. Theys, T., Fati, D., Schrevens, O., Tarnowska, A., Zienkiewicz, M., 2014. From lab to
Hemalatha, K., Rani, K.Y., 2017. Multiobjective optimization of unseeded and seeded plant: first industrial experience of the new high efficiency dihydrate hemihy-
batch cooling crystallization processes. Ind. Eng. Chem. Res. 56 (20), 6012–6021. drate process (DA-HF) process. Procedia Eng. 83, 181–187. doi:10.1016/[Link].
doi:10.1021/[Link].7b00586. 2014.09.036.
Holmes, H., Mesmer, R., 1999. Isopiestic studies of H3PO4 (aq) at elevated tempera- Wang, W., Zeng, D., Chen, Q., Yin, X., 2013. Experimental determination and
tures. J. Solution Chem. 28 (4), 327–340. doi:10.1023/A:1022651710922. modeling of gypsum and insoluble anhydrite solubility in the system
Irish, D., Chen, H., 1971. Raman spectral study of bisulfate-sulfate systems. ii. consti- CaSO4–H2SO4–H2O. Chem. Eng. Sci. 101, 120–129. doi:10.1016/[Link].2013.06.
tution, equilibriums, and ultrafast proton transfer in sulfuric acid. J. Phys. Chem. 023.
75 (17), 2672–2681. doi:10.1021/j100686a024. Wu, S., McLean, K.A., Harris, T.J., McAuley, K.B., 2011. Selection of optimal pa-
Joao, S., Durand, A., Schrevens, O., 2016. Plant operability optimization through dy- rameter set using estimability analysis and MSE-based model-selection crite-
namic simulation, a case study focused on phosphoric acid concentration unit. rion. Int. J. Adv. [Link]. 3 (3), 188–197. doi:10.1504/IJAMECHS.2011.
Procedia Eng. 138, 378–389. doi:10.1016/[Link].2016.02.097. 042615.

14
I. Bouchkira, A.M. Latifi, L. Khamar et al. Computers and Chemical Engineering 156 (2022) 107536

Yang, H., Zhao, Z., Zeng, D., Yin, R., 2016. Isopiestic measurements of water activity Yao, K.Z., Shaw, B.M., Kou, B., McAuley, K.B., Bacon, D., 2003. Modeling ethy-
for the H2SO4–H3PO4–H2O system at 298.15 K. J. Solution Chem. 45 (11), 1580– lene/butene copolymerization with multi-site catalysts: parameter estimabil-
1587. doi:10.1007/s10953- 016- 0516- 4. ity and experimental design. Polym. React. Eng. 11 (3), 563–588. doi:10.1081/
Yang, L., Cao, J., Luo, T., 2018. Effect of Mg2+, Al3+, and Fe3+ ions on crystallization PRE-120024426.
of type α hemi-hydrated calcium sulfate under simulated conditions of hemi- Zdanovskii, A., Vlasov, G., 1968. Determination of the boundaries of the reciprocal
hydrate process of phosphoric acid. J. Cryst. Growth 486, 30–37. doi:10.1016/j. transformation of CaSO4·2H2O and γ -CaSO4 in H2SO4 solutions. Russ. J. Inorg.
jcrysgro.2018.01.014. Chem. 13, 1318–1319.

15