# Multivariate optimization and inverse design of multiphase reactive systems via deep learning and 3D CFD-PBM simulations

Antonello Raponi, Daniele Marchisio

► To cite this version:

Antonello Raponi, Daniele Marchisio. Multivariate optimization and inverse design of multiphase reactive systems via deep learning and 3D CFD-PBM simulations. Chemical Engineering Journal, 2026, 531, pp.173908. 10.1016/j.cej.2026.173908. emse-05523779

HAL Id: emse-05523779

https://hal-emse.ccsd.cnrs.fr/emse-05523779v1

Submitted on 23 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Multivariate optimization and inverse design of multiphase reactive systems via deep learning and 3D CFD-PBM simulations

Antonello Raponi<sup>a,b,10,*</sup>, Daniele Marchisio<sup>b</sup>

## ARTICLE INFO

**Keywords:**

- Multivariate optimization
- Artificial Intelligence
- Population balance model
- Computational fluid dynamics
- Digital twins
- Inverse design

## ABSTRACT

This work presents a deep learning-driven framework for multivariate optimization, focusing on the identification of precipitation kinetics parameters in a complex three-dimensional (3D) Computational Fluid Dynamics-Population Balance Model (CFD-PBM) model describing Mg(OH)<sub>2</sub> precipitation in T- and Y-mixer static reactors. A numerical dataset was first generated using the CFD-PBM model to train an inverse design deep learning neural network (mirror model), which takes characteristic particle dimensions (\dagger) as input and predicts the corresponding kinetic parameters (\phi) as output. Experimental particle size distributions (PSD) from the T-mixer were then provided to the trained mirror model to predict the kinetic parameter vector. The predicted parameters were subsequently validated by comparing the CFD-PBM outputs with experimental PSD from the Y-mixer configuration. The framework demonstrates excellent generalization, accurately reproducing PSD for conditions not included in the training set. Moreover, the same numerical dataset was also employed to train a surrogate CFD-PBM model which takes kinetic parameters (\phi) as input and returns characteristic particle dimensions (\dagger) as output. The surrogate CFD-PBM model was subsequently used to perform a sensitivity analysis of the kinetic parameters. Minor deviations are observed in the estimation of the smallest (d<sub>10</sub>) and largest (d<sub>43</sub>) particle size fractions, consistent with physical aggregation and collision mechanisms. A key outcome of this study is that the proposed framework enables, for the first time, the use of fully 3D CFD-PBM simulations equipped with a rigorously calibrated kinetic model, without requiring empirical simplifications or iterative fitting procedures. This significantly enhances the predictive capability of 3D CFD-PBM models in nonlinear reactive systems, providing a computationally efficient and generalizable tool for kinetic parameter identification in precipitation and crystallization processes.

## 1. Introduction

Magnesium has recently been classified among the Critical Raw Materials (CRMs), highlighting the importance of efficient recovery and valorization strategies [1,2]. One effective approach involves using magnesium to produce magnesium hydroxide, Mg(OH)<sub>2</sub>, a versatile compound with applications spanning flame-retardant fillers, acidic waste neutralization, pharmaceuticals, and pulp and paper preservation [3–9]. Traditionally, Mg(OH)<sub>2</sub> has been synthesized through hydrothermal processes [10], but the high energy requirements have motivated a shift towards precipitation-based methods, which are more sustainable and cost-effective [11].

Precipitation involves the reaction of magnesium sources such as brines or bitterns with alkaline solutions, leading to supersaturation and particle formation when the solubility limit is exceeded [12]. The resulting Particle Size Distribution (PSD) is influenced by primary

nucleation, molecular growth, and aggregation phenomena. Population Balance Models (PBM) provide a powerful framework to track PSD evolution and to optimize process conditions for desired particle characteristics [13,14]. Accurate determination of the kinetic parameters governing these processes is therefore essential for process optimization [15,16].

Traditionally, model parameters are identified by fitting model predictions to experimental data using optimization algorithms. Classical methods, such as the Conjugate Gradient (CG) method, are effective when the underlying first-principles model can be evaluated with relatively low computational cost. Similarly, heuristic and meta-heuristic approaches, including Particle Swarm Optimization (PSO) [17] or Crow Search Optimization (CSO) [18], provide robust alternatives for exploring high-dimensional search spaces [19–22]. However, these methods can become impractical when the first-principle model is highly

<footnote><sup>a</sup> Centre SPIN, CNRS UMR 5307 Laboratoire Georges Friedel, Mines Saint-Etienne, France</footnote>

<footnote><sup>b</sup> Politecnico di Torino, Institute of Chemical Engineering, Department of Applied Science and Technology, 10129 Torino, Italy</footnote>

<footnote>* Corresponding author at: Centre SPIN, CNRS UMR 5307 Laboratoire Georges Friedel, Mines Saint-Etienne, France.
E-mail address: antonello.raponi@emse.fr (A. Raponi).</footnote>
non-linear, three-dimensional (3D), and computationally expensive, as the numerous repeated evaluations required for optimization may exceed reasonable computing time, even on high-performance hardware. Surrogate-Based Optimization (SBO) has emerged as a strategy to alleviate the computational burden in such cases. SBO replaces the original high-fidelity model with a computationally cheaper surrogate, constructed from a low-fidelity approximation combined with a correction procedure to increase accuracy. This approach has been successfully applied in engineering contexts and in marine ecosystem models of NPZD type (where N stands for nutrients, P for phytoplankton, Z for Zooplankton and D for Detritus). In particular, it has been extended from one-dimensional applications to a 3D biogeochemical model, where the low-fidelity component is represented by a simulation with a reduced number of spin-up iterations, and a multiplicative correction operator is used to extrapolate its output towards the high-fidelity solution [23]. While SBO significantly reduces computational time and can accurately estimate most parameters, it does not overcome the inherent limitations of traditional optimization methods, such as sensitivity to the initial guess, limited exploration of multi-objective landscapes, or reliance on the quality of available data.

These limitations motivate the use of alternative strategies that combine first-principles models with data-driven approaches, such as deep learning-assisted inverse design [24–26]. By leveraging neural networks trained on numerically generated model outputs, it becomes possible to map observed system outputs (e.g., PSD) directly to the underlying kinetic parameters, achieving both computational efficiency and improved robustness in parameter identification. These methods enable efficient exploration of parameter spaces in complex systems, providing accurate predictions with reduced computational cost. Recent studies have demonstrated that machine learning can serve as an efficient surrogate for computationally expensive CFD simulations. Comparisons between neural networks, tree-based models (Random Forest, XGBoost), and other regressors show that alternative algorithms can achieve comparable predictive accuracy while reducing computational cost [27–29]. In addition, recent comprehensive reviews have highlighted the growing role of machine learning in multiphase flows and reactor modeling. For example, Zhu et al. [30] provide an extensive survey of machine learning applications to hydrodynamics, transport phenomena, and reaction kinetics in multiphase flows and reactors, discussing both physics-based and data-driven approaches and outlining opportunities and challenges in this domain.

In this work, we employ a fully three-dimensional Computational Fluid Dynamics–Population Balance Model (CFD–PBM) to describe the evolution of Mg(OH)₂ particle sizes, accounting for homogeneous and heterogeneous nucleation, molecular growth, and irreversible aggregation via cementation. The PBM is solved using the Quadrature Method of Moments (QMOM) [31,32], and the resulting high-fidelity numerical data are used to train a fully connected deep learning neural network for inverse design. Experimental measurements guide the deep learning neural network in identifying the unknown kinetic parameters directly from the characteristic particle dimensions, avoiding empirical simplifications and eliminating the need for traditional iterative fitting procedures. This approach enables, for the first time, the development of CFD–PBM simulations equipped with a rigorously calibrated kinetic model for a complex reactive, multiphase, and strongly nonlinear system. In addition, it enables a framework with potential for digital twin applications, capable of predictive simulations and real-time response through regression neural networks. Two experimental datasets from our previous studies [33,34] are employed: one for tuning the model and one for independent validation, allowing to assess the generalization capability of the proposed framework. The manuscript is organized as follows: Section 2 describes the experimental protocols and datasets for the T- and Y-mixer reactors; Section 3 details the PBM and its key inputs and outputs; Section 4 introduces the deep learning-driven parameter identification methodology; Section 5 presents the results and discusses the performance of the inverse-design strategy; and Section 6 summarizes the main findings and outlines future extensions.

## 2. Methodology

This section describes the experimental setups, the main steps for acquiring experimental data and the operating conditions.

**Experimental setup.** Some of our previous publications [33,34] describe a well-tuned, and validated experimental procedure for stabilizing and characterizing Mg(OH)₂ suspensions. This procedure was applied to suspensions obtained through the two experimental apparatuses, T- and Y- static mixers, shown in Fig. 1.

The Y-mixer apparatus is composed by a Y-mixer, followed by two diverging channels and a final pipe of constant diameter. Although the procedure details have been widely illustrated in the aforementioned publications [33,34], they are here briefly reviewed for clarity. Magnesium chloride (MgCl₂) solution is fed in one of the mixer's arms and the sodium hydroxide (NaOH) solution in the other. After the impact in the mixing channel, the two solutions react to form a Mg(OH)₂ suspension according to the following reaction:

$$ \mathrm{MgCl}_{2}(\mathrm{aq}) + 2\mathrm{NaOH}(\mathrm{aq}) \rightarrow \mathrm{Mg(OH)}_{2}(\mathrm{s}) \downarrow \quad (1) $$

The collected Mg(OH)₂ suspension was stabilized and then analyzed by Dynamic Light Scattering (DLS) measurements to obtain the PSD. The investigated operating conditions are reported in Table 1:

In Table 1 only the Mg²⁺ is reported concentration is reported, as the feed streams respect the stoichiometric ratio of 2:1. The characteristic dimensions both for the T- and Y-mixers are presented in Fig. 2.

## 3. Population balance model

In this section, we present the CFD-PBM model used to generate the numerical dataset for training the deep learning neural network. This model assumes a one-way coupling approach. It means that the particles have no inertia and follow the flow field. This hypothesis is justified by the dilute operating conditions investigated, for which the solid volume fraction remains low and the suspension properties (density and viscosity) are close to those of the liquid phase. As a consequence, the feedback of particle formation and growth on the velocity and turbulence fields can be reasonably neglected. Moreover, the characteristic timescales associated with nucleation and particle growth are significantly larger than those governing the establishment of the flow field, supporting the decoupling between hydrodynamics and particle population dynamics. Although the flowchart has been introduced in our previous publication [35], we briefly present it here. First, the CFD simulation is run to solve the flow (ū) and turbulence (k-ε) fields. A second-order discretization scheme was applied to both the convective and diffusive terms in the momentum and turbulence transport equations to ensure higher numerical accuracy and reduce numerical diffusion. Fig. 3 shows the meshes wireframes used in the simulations. For the sake of visualization both meshes were cropped and enlarged:

In Table 2 are reported the boundary conditions used in the simulations:

Once these fields are known, the equations for the macro- and micro-mixing are solved. In this regard, the macro-mixing is represented by the evolution of a non-reactive scalar, namely the mixture fraction (ā), indicating the interaction of the fluids at the macro-scale. The governing equation is as follows:

$$ \frac{\partial \bar{\alpha}}{\partial t} + \bar{u} \cdot \nabla \bar{\alpha} = D_t \nabla^2 \bar{\alpha} \quad (2) $$

where D_t is the turbulent diffusion. The micro-mixing is represented by the evolution of the variance of the mixture fraction (ā'²), indicating the interaction of the fluids at the molecular level. The transport equation for the variance of the mixture fraction is:

$$ \frac{\partial \bar{\alpha}'}{\partial t} + \bar{u} \cdot \nabla \bar{\alpha}^{\prime 2} = D_t \nabla^2 \bar{\alpha}^{\prime 2} + 2 D_t (\nabla \bar{\alpha} \cdot \nabla \bar{\alpha}) - \frac{C_{\phi} \epsilon}{2 k} \bar{\alpha}^{\prime 2} \quad (3) $$
<figure><img src="image_3.png" /><figcaption>Fig. 1. Experimental setups: T-mixer (top) and Y-mixer (bottom).</figcaption></figure>

<table><caption>**Table 1**

Operating conditions for the T- and Y-mixer.</caption><thead><tr><th rowspan="2">Dataset</th><th rowspan="2">Geometry</th><th colspan="5">Operating conditions</th></tr><tr><th>Concentration (M)</th><th>Residence time (ms)</th><th>Reynolds number (-10<sup>4</sup>)</th><th>Flow rate (mL/min)</th><th>Estimate mixing time (ms)</th></tr></thead><tbody><tr><td rowspan="5">#1</td><td rowspan="5">T-Mixer</td><td>0.125</td><td>~3</td><td>2.7</td><td>2320</td><td>2</td></tr><tr><td>0.25</td><td>~3</td><td>2.7</td><td>2320</td><td>2</td></tr><tr><td>0.5</td><td>~3</td><td>2.7</td><td>2320</td><td>2</td></tr><tr><td>0.75</td><td>~3</td><td>2.7</td><td>2320</td><td>2</td></tr><tr><td>1</td><td>~3</td><td>2.7</td><td>2320</td><td>2</td></tr><tr><td rowspan="6">#3</td><td rowspan="6">Y-Mixer</td><td>0.01</td><td>~9.10<sup>3</sup></td><td>1.7</td><td>835</td><td>0.6</td></tr><tr><td>0.025</td><td>~9.10<sup>3</sup></td><td>1.7</td><td>835</td><td>0.6</td></tr><tr><td>0.05</td><td>~9.10<sup>3</sup></td><td>1.7</td><td>835</td><td>0.6</td></tr><tr><td>0.125</td><td>~3.6·10<sup>2</sup></td><td>1.7</td><td>835</td><td>0.6</td></tr><tr><td>0.5</td><td>~3.6·10<sup>2</sup></td><td>1.7</td><td>835</td><td>0.6</td></tr><tr><td>1</td><td>~3.6·10<sup>2</sup></td><td>1.7</td><td>835</td><td>0.6</td></tr></tbody></table>

<table><caption>**Table 2**

Boundary conditions used for the CFD simulations in the T- and Y-mixers.</caption><thead><tr><th>Boundary</th><th>Velocity</th><th>Pressure</th><th>Turbulence fields</th></tr></thead><tbody><tr><td>Inlet 1</td><td>Specified flow rate</td><td>Neumann</td><td>Correlation</td></tr><tr><td>Inlet 2</td><td>Specified flow rate</td><td>Neumann</td><td>Correlation</td></tr><tr><td>Outlet</td><td>Neumann</td><td>0</td><td>Neumann</td></tr><tr><td>Walls</td><td>No-slip</td><td>Neumann</td><td>Wall function</td></tr></tbody></table>

In addition to the accumulation, advection, and diffusion terms, Eq. (3) contains two other contributions that merit closer analysis. The first is the production term, $2D_t(\nabla\bar{\alpha} \cdot \nabla\bar{\alpha})$, which represents the macroscopic fluid interaction responsible for generating variance. The second is the dissipation term associated with the smallest scales (i.e., the Batchelor scale), where the variance is dissipated. Unlike the production term, the dissipation term cannot be directly determined and thus requires a closure relationship [36]. In our previous studies [33,34], the parameter $C_\phi$ was assigned a value of 2. The mixture fraction and its variance enable the model to account for the dampening effect of mixing on the chemical reaction. Specifically, $Mg^{2+}$ and $OH^-$ must undergo both macro- and micro-mixing before reacting according to Eq. (1), which slows down the precipitation process. In Raponi et al. [33] (Supplementary Information), we derived an analytical expression linking the concentration of reactants entering the system to the fraction that is actually able to react at the microscale to form $Mg(OH)_2$. For clarity, only the final form of this equation is reported here:

$$[Mg(OH)_2]_{(aq)} = \min \left( \left[ Mg^{2+} \right], \frac{[OH^-]}{2} \right) f \left( \bar{\alpha}, \bar{\alpha}'^2, \bar{\alpha}_s \right) \quad (4)$$

where $Mg^{2+}$ and $OH^-$ represent the concentrations of the respective species in the system, and $f(\bar{\alpha}, \bar{\alpha}'^2, \bar{\alpha}_s)$ is a weighting function that accounts for the effect of mixing. Here, $\bar{\alpha}_s$ denotes the stoichiometric mixture fraction, defined as:

$$\bar{\alpha}_s = \frac{2 [Mg^{2+}]^{(in)}}{2 [Mg^{2+}]^{(in)} + [OH^-]^{(in)}} \quad (5)$$

where the superscript “(in)” denotes the reactant concentration at the inlets. In the mixture fraction framework, it corresponds to the mixture fraction at which the reactants are present in stoichiometric proportions. By assuming that the $Mg(OH)_2_{(aq)}$ in the aqueous phase precipitates instantaneously as a solid, the supersaturation can be defined as:

$$S = \frac{\gamma_{\pm}^3 [Mg^{2+}]_{(aq)} [OH^-]_{(aq)}^2}{k_{sp}} - 1 \quad (6)$$

where $[Mg^{2+}]_{(aq)}$ and $[OH^-]_{(aq)}$ are the concentrations calculated from Eq. (4) according to stoichiometry, $\gamma_{\pm}$ is the activity coefficient, and the solubility product: $k_{sp} = 10^{-10.88}$. The build-up of supersaturation thus initiates the precipitation process. In this work, we consider primary nucleation, both homogeneous and heterogeneous, denoted as J:

$$J = A_1 e^{\left( -\frac{B_1}{\ln^2(S+1)} \right)} + A_2 e^{\left( -\frac{B_2}{\ln^2(S+1)} \right)} \quad (7)$$

molecular growth, G:

$$G = k_g S^g \quad (8)$$
<figure><img src="image_4.png" /><figcaption>Fig. 2. Visualization of the experimental characteristic dimensions.</figcaption></figure>

and the aggregation rate, $\beta_{agg}$:

$$ \beta_{agg} = \beta_{col}\eta_{agg} \qquad (9) $$

The aggregation rate is proportional to the collision frequency, $\beta_{col}$, which represents the number of collisions per unit time between primary particles, and to the sticking probability, $\eta_{agg}$, describing the likelihood that a collision results in actual aggregation (i.e., whether the primary particles bind together or not). The collision frequency incorporates contributions from both Brownian motion and turbulence.

$$ \beta_{col} = 10^{C_1} (\beta_{br} + \beta_{tr}) \qquad (10) $$

The contribution due to Brownian motion is described by the following equation:

$$ \beta_{br} = \frac{2k_B T (L + \lambda)^2}{3\mu L\lambda} \qquad (11) $$

In this expression, $k_B$ denotes the Boltzmann constant, $T$ is the fluid temperature, $\mu$ represents the dynamic viscosity of the fluid, and $L$ and $\lambda$ are the sizes of the colliding particles. The turbulent contribution, in contrast, is given by:

$$ \beta_{tr} = \sqrt{\frac{8\pi}{15} \sqrt{\frac{\epsilon}{\nu} (L + \lambda)^3}} \qquad (12) $$

The sticking probability, on the other hand, can be modeled as follows:

$$ \eta_{agg} = e^{-\theta} \qquad (13) $$

Here, $\theta$ denotes the ratio between the cementation time, $t_{cem}$, and the interaction time, $t_{int}$. The interaction time is assumed to correspond to the characteristic timescale of the Kolmogorov micro-scale:

$$ t_{int} = \sqrt{\frac{\nu}{\epsilon}} \qquad (14) $$

and corresponds to the duration during which the two particles remain close enough to interact, while the cementation time is defined as:

$$ t_{cem} = \frac{D_b}{s(\delta)G} \qquad (15) $$

$$ D_b = \frac{L_{eqv} \rho_p^{0.5} (\epsilon\nu)^{0.25}}{A_p^{0.5}} \qquad (16) $$

$$ L_{eqv} = \frac{L\lambda}{(L^2 + \lambda^2 - L\lambda)^{0.5}} \qquad (17) $$

It corresponds to the time needed for the irreversible cementation of two primary particles. The function $s(\delta)$ is a shape factor, defined as follows [37]:

$$ s(\delta) = \frac{4(1 + \delta - \delta')}{1/3 + \delta - \delta' - (\delta - \delta')^2 (2\delta/3 + \delta'/3)} \qquad (18) $$

$$ \delta = L/\lambda \qquad (19) $$

$$ \delta' = \sqrt{\delta^2 - 1} \qquad (20) $$

Eqs. (7), (8), (10), and (16) involve eight tuning parameters ($A_1$, $B_1$, $A_2$, $B_2$, $k_g$, $g$, $C_1$, and $A_p$). The parameters $A_i$ and $B_i$ control the intensity of primary nucleation, $k_g$ regulates the rate of molecular growth, and $g$ characterizes the growth mechanism (e.g., $g = 1$ indicates diffusion-controlled growth). The parameter $C_1$ modulates the collision frequency, while $A_p$ is proportional to the strength of the crystalline bridge between primary particles. The evolution of primary particles resulting from the precipitation phenomena described above was tracked using the Quadrature Method of Moments (QMOM) [38,39]. The transport equation for the $k$th moment is given by:

$$ \frac{\partial \bar{m}_k}{\partial t} + \bar{u} \cdot \nabla \bar{m}_k = D_t \nabla^2 \bar{m}_k + J L_c^k + k G m_{k-1} + \bar{B}_k - \bar{D}_k \qquad (21) $$

where $L_c$ denotes the critical size of a stable nucleus (1 nm), and $\bar{B}_k$ and $\bar{D}_k$ represent the birth and death terms due to particle aggregation. In the present study, three quadrature nodes (corresponding to six moments) were employed. Finally, the population balance model (PBM) is closed by the mass balances of the reacting species, $[Mg^{2+}]$ and $[OH^{-}]$:

$$ \frac{\partial \bar{c}_i}{\partial t} + \bar{u} \cdot \nabla \bar{c}_i = D_t \nabla^2 \bar{c}_i - \frac{\nu_i \rho_p k_v}{M_p} (J L_c^3 + 3 G m_2) \qquad (22) $$

where $\nu_i$ denotes the stoichiometric coefficient (i.e., $\nu_{Mg^{2+}} = 1$ and $\nu_{OH^{-}} = 2$), $\rho_p$ is the particle density, $M_p$ is the molecular weight and $k_v = \pi/6$ is the particle shape factor.

## 4. Parameters identification

In this section, we present the application of a deep learning-aided optimization framework to identify the kinetic parameters of a complex, highly non-linear, 3D CFD-PBM model. Fig. 4 helps to visualize the dataset generation and the inverse design (mirror model) used.
<figure><img src="image_5.png" /><figcaption>Fig. 3. T- and Y-mixer mesh wireframe.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Dataset generation and inverse design framework.</figcaption></figure>

The vector $\vec{\phi}$ represents the kinetic parameters ($\vec{\phi} = \phi(A_1, A_2, B_1, B_2, k_g, g, C_1, A_p)$), while the vector $\vec{d}$ denotes the characteristic dimensions ($\vec{d} = d(d_{10}, d_{21}, d_{32}, d_{43})$) where the characteristic dimension $d_{ij}$ is defined as the ratio between the outlet mean value of the $i$th order moment and the outlet mean value of the $j$th order moment [13]:

$$d_{ij} = \left( \frac{\vec{m}_i}{\vec{m}_j} \right)_{\text{outlet}} \qquad (23)$$

The first step involves generating a numerical dataset using the CFD-PBM model (Fig. 4 left). Specifically, the CFD-PBM model takes as input the kinetic parameters and, under given operating conditions, returns the corresponding characteristic dimensions. This dataset is designed to remain consistent with the physics of the problem. The kinetic parameters are randomly sampled within their physically meaningful domain, as readily available in the literature [33,40]. The lower and upper bounds for each kinetic parameter, consistent with our previous publications and standard handbooks, are reported in Table 3.

To provide a sufficiently rich dataset for training the mirror model, 400 CFD-PBM simulations were generated by sampling the kinetic parameters across their physically meaningful ranges. Each simulation was performed under one of the five experimental $Mg^{2+}$ concentrations in the T-mixer, meaning that the same operating condition was simulated multiple times with different kinetic parameter sets. This procedure allows the CFD-PBM model to explore a wide range of physically

<table><caption>Table 3
Parameter constraints and units.</caption><thead><tr><th></th><th>A<sub>1</sub></th><th>A<sub>2</sub></th><th>B<sub>1</sub></th><th>B<sub>2</sub></th><th>k<sub>g</sub></th><th>g</th><th>C<sub>1</sub></th><th>A<sub>p</sub></th></tr></thead><tbody><tr><td>Lower bound</td><td>10<sup>19</sup></td><td>10<sup>10</sup></td><td>200</td><td>10</td><td>10<sup>-15</sup></td><td>1</td><td>0</td><td>10<sup>0</sup></td></tr><tr><td>Upper bound</td><td>10<sup>29</sup></td><td>10<sup>18</sup></td><td>400</td><td>10<sup>2</sup></td><td>10<sup>-9</sup></td><td>2</td><td>1</td><td>10<sup>7</sup></td></tr><tr><td>Units</td><td colspan="2">particle no.</td><td colspan="2">m</td><td colspan="2">m</td><td colspan="2">N</td></tr><tr><td></td><td colspan="2">m<sup>3</sup>s<sup>1</sup></td><td colspan="2">s</td><td colspan="2">s</td><td colspan="2">m<sup>2</sup></td></tr></tbody></table>

plausible system behaviors and effectively cover the kinetic parameter space. The numerical dataset generated by the CFD-PBM model is subsequently used to train a deep learning neural network (Fig. 4 right). The approach is referred to as inverse design (mirror model) because the deep learning neural network receives as input the vector of characteristic dimensions ($\vec{d}$) and outputs the corresponding kinetic parameters ($\vec{\phi}$). Starting from the numerical dataset, the deep learning neural network is trained and tested. Following the training and testing phases, the experimental characteristic dimensions obtained from the T-mixer are provided to the deep learning neural network, which in turn predicts the kinetic parameters that, according to the first-principles model used for the dataset generation, best reproduce the experimental characteristic dimensions. This procedure is repeated for all experimental conditions (five concentrations in the case of the T-mixer). Since each experimental set of characteristic dimensions yields a predicted
<table><caption>**Table 4**
Summary of the optimized deep learning neural network architectures obtained from the Keras Tuner for the two CFD-PBM datasets.</caption><thead><tr><td>Dataset</td><td>Hidden layers</td><td>Neurons per layer</td><td>Learning rate</td></tr></thead><tbody><tr><td>200 simulations</td><td>3</td><td>192-68-112</td><td>1e-3</td></tr><tr><td>400 simulations</td><td>3</td><td>192-160-160</td><td>1e-3</td></tr></tbody></table>

parameter set, the kinetic parameters are subsequently averaged ($\bar{\phi}$), and their standard deviation is calculated ($\bar{\phi}'^2$). It is worth noting that this averaging procedure is not intended to provide a statistically converged estimate in the classical sense. Rather, each experimental operating condition yields an independent solution, corresponding to a kinetic parameter set that reproduces the measured characteristic dimensions. The mean kinetic parameter vector ($\bar{\phi}$) is representative of a set that best captures the ensemble of experimental conditions, while the standard deviation ($\bar{\phi}'^2$) provides a quantitative measure of the robustness and consistency of the parameter identification procedure with respect to changes in operating conditions. It is important to emphasize that the numerical dataset generated from the CFD-PBM model is used solely for training the inverse model and does not aim to reproduce the specific experimental datasets. The broad sampling of the kinetic parameter space is required to ensure that the deep learning model learns a well-posed inverse relationship between characteristic dimensions and kinetic parameters, rather than overfitting a narrow region of the parameter space. The experimental data are introduced only at the inference stage of the trained mirror model. Each experimental operating condition is treated as an independent inverse problem, yielding a kinetic parameter set that reproduces the corresponding experimental characteristic dimensions based on the first-principles model. The averaging of the identified parameters across experimental conditions therefore provides a representative kinetic parameter vector that consistently captures the ensemble of experiments, while the associated standard deviation quantifies the robustness of the identification with respect to changes in operating conditions. At this stage, two aspects are crucial: (i) the choice of the deep learning neural network architecture, and (ii) the extension of the numerical dataset, i.e., the number of simulations required for the mirror model to adequately learn the correlations between the characteristic-dimension vector and the kinetic-parameter vector. To select the optimal architecture and training parameters, we employed the Keras Tuner, which systematically explores hyperparameters including the number of hidden layers, neurons per layer, and learning rate. This ensures a reproducible and data-driven selection of hyperparameters, rather than arbitrary manual choices. Two numerical datasets were used for training the mirror model. The first dataset consisted of 200 CFD-PBM simulations, while the second included 400 CFD-PBM simulations. In both cases, the data were split into 90% for training and 10% for testing. Following the same approach adopted in our previous work [24], the Keras Tuner was executed for both datasets. The outcomes of the two optimization processes are summarized in Table 4.

The analysis showed that the model trained with 200 CFD-PBM simulations exhibited large variability in the predicted parameters, indicating insufficient generalization. Conversely, the model trained with 400 simulations provided stable and consistent predictions, with significantly reduced standard deviations.

## 5. Results and discussion

In this section, we present the results of this multivariate optimization procedure. Once the mean values ($\bar{\phi}$) and standard deviations ($\bar{\phi}'^2$) of the kinetic parameters were obtained, two parallel analyses were performed: the former aimed at validating the overall framework, and the latter focused on estimating the intrinsic uncertainty of the kinetic parameter identification procedure. The validation of the numerical

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Comparison between the experimental data (red dashed line), surrogate model (gray dots) and CFD-PBM predictions (black dashed line).</figcaption></figure>

framework was carried out by simulating the operating conditions reported in Fig. 2, using the mean kinetic parameters predicted by the inverse design (mirror model) both for the T- and Y-mixer configurations. Subsequently, a surrogate CFD-PBM model was trained using the dataset of 400 simulations for the T-mixer. This surrogate model, which maps the kinetic parameter vector to the corresponding characteristic dimensions, was employed to assess both the robustness of the framework and the uncertainty associated with the predicted kinetic parameters. Specifically, the surrogate CFD-PBM model was used to predict the characteristic dimensions for 1000 randomly sampled kinetic parameter vectors, drawn within the range of the mean value ± the standard deviation ($\bar{\phi} - \bar{\phi}'^2$ and $\bar{\phi} + \bar{\phi}'^2$). The mean and standard deviation of the predicted characteristic dimensions were then computed for all the operating conditions. Both the CFD-PBM predictions and the uncertainty analysis are reported in Fig. 5.

The CFD-PBM model, applied to both the T- and Y-mixer configurations, employed the mean kinetic parameter vector $\bar{\phi}$ predicted by the
<figure><img src="image_8.png" /><figcaption>Fig. 6. Contour plot of the $d_{10}$ in the Y-mixer at three $Mg^{2+}$ concentrations: 0.01 M (left), 0.05 M (center), 1 M (right).</figcaption></figure>

inverse design. As shown in Fig. 5, the use of these mean kinetic parameters allows the CFD-PBM to accurately capture the overall system dynamics, demonstrating a high level of agreement with the experimental data. The sensitivity analysis for the kinetic parameters further confirms the reliability of the framework. In particular, it successfully reproduces the trends observed in the CFD-PBM simulations for the T-mixer, even when evaluated on kinetic parameters set that were not included in the training phase, namely, those obtained from the inverse design (i.e., 1000 random sets included between $\overline{\overline{\phi}} - \overline{\overline{\phi}}'^2$ and $\overline{\overline{\phi}} + \overline{\overline{\phi}}'^2$). This highlights the strong generalization capability of the surrogate CFD-PBM model. It should be noted that the surrogate CFD-PBM model was trained exclusively using the dataset generated for the T-mixer configuration. For the Y-mixer, the kinetic parameters identified from the T-mixer were directly applied to run the one-way coupled CFD-PBM simulations. Since only six CFD-PBM simulations were available for the Y-mixer, corresponding to the experimental operating points, the construction of a surrogate model capable of reliably mapping a wide range of kinetic parameters to characteristic dimensions was not feasible. Consequently, the validation for the Y-mixer relies solely on the one-way coupled CFD-PBM simulations performed using the identified mean kinetic parameters. From Fig. 5, it can also be observed that the kinetic parameters estimated through the inverse design tend to slightly underestimate the characteristic dimension $d_{10}$ for both the surrogate and the full CFD-PBM model. Moreover, the surrogate model exhibits slightly wider confidence intervals for the higher concentration cases, particularly for the characteristic dimension $d_{43}$. This trend can be intuitively explained from a physical standpoint. The surrogate model is inherently derived from the CFD-PBM dataset used for training, and the CFD-PBM model itself predicts a higher number of particles at increasing concentrations. The $d_{43}$ metric represents the right tail of the PSD, corresponding to the larger particles. Therefore, at higher concentrations, the increased number of particle collisions and the formation of larger aggregates lead to a broader distribution and, consequently, greater variability in $d_{43}$. Conversely, the standard deviation of $d_{10}$ remains nearly constant with concentration, indicating that the formation of smaller particles is less sensitive to concentration-driven effects. To further elucidate the influence of local mixing on particle formation, Fig. 6 reports the contour plots of the characteristic dimensions $d_{10}$ within the Y-mixer for three selected $Mg^{2+}$ concentrations.

These concentrations were chosen based on the trends observed in Fig. 5, namely a low concentration (0.01 M), a value close to the numerical minimum (0.05 M), and the highest concentration considered

(1 M). The inclusion of these contour plots is of practical relevance, since the proposed inverse-design framework enables the optimization of kinetic parameters even for highly nonlinear and fully 3D systems. Consequently, it becomes possible to examine the internal behavior of the reactor and gain insight into how spatial gradients in mixing affect the early stages of the precipitation process. All simulations were carried out under identical hydrodynamic conditions, with only the $Mg^{2+}$ concentration being varied. This ensures that the mixing characteristics remain unchanged across the cases, as both $\bar{a}$ and $\bar{a}'^2$ were kept constant, together with a stoichiometric mixing ratio, $\bar{a}_s$, equal to 0.5. As shown by Eq. (5), the driving force increases with the concentration of the reactive species. The contour plots of $d_{10}$ clearly highlight how this effect interacts with the spatially non-uniform mixing field inside the Y-mixer. At the lowest concentration (0.01 M), the driving force is insufficient to trigger nucleation within the mixing section. As a result, nucleation mainly occurs in the diverging channels and in the downstream tube, where particles experience more homogeneous turbulent conditions and longer residence times. This leads to relatively uniform and larger final particle sizes. A different behavior is observed at the intermediate concentration (0.05 M). Here, the increased driving force promotes earlier nucleation, which now takes place within the Y-mixer itself, as visible in Fig. 6 (center). In this region, the strong local turbulence enhances the formation of a larger number of nuclei, accelerating the local consumption of supersaturation in a non-uniform manner. Consequently, less supersaturation remains available for subsequent aggregation, resulting in smaller particles compared to the regime where nucleation occurs downstream. At the highest concentration (1 M), the mechanisms observed at 0.05 M become more pronounced. While nucleation still begins within the Y-mixer, the significantly larger reactant concentration ensures a higher residual supersaturation at equivalent spatial coordinates, which favors aggregation. This yields a qualitative trend similar to the intermediate case but leads to quantitatively larger particles, consistent with the stronger driving force. Overall, these contour plots underline the critical interplay between spatial mixing heterogeneities and local supersaturation in determining the onset of particle formation. More importantly, they demonstrate the capability of the proposed inverse-design-based 3D CFD-PBM framework to generate spatially resolved and physically interpretable predictions. This establishes a solid foundation for future digital twins able to more accurately represent real systems and to support practical applications such as process control and reactor design.
# 6. Conclusions

In this study, a deep learning (DL)-driven framework was successfully developed and applied for the identification of kinetic parameters in a complex CFD-PBM model describing $Mg(OH)_2$ precipitation in static mixer reactors. The proposed methodology combines first-principles modeling with data-driven optimization, leveraging an inverse design (mirror model) deep learning neural network trained on a numerically generated CFD-PBM dataset. The inverse design approach enabled the identification of the kinetic parameter vector directly from experimentally measured characteristic particle sizes, significantly reducing the computational cost typically associated with iterative parameter fitting. The mirror model trained on 400 CFD-PBM simulations demonstrated superior generalization capability and parameter stability compared to the smaller dataset of 200 simulations, confirming the importance of dataset extension in capturing the nonlinear correlations between kinetics and morphology. The mean kinetic parameters predicted by the mirror model accurately reproduced the experimental trends for both the T- and Y-mixer configurations, validating the robustness of the coupled CFD-PBM-DL framework. Sensitivity and uncertainty analyses further revealed that the surrogate model reliably captures the system dynamics even for unseen kinetic configurations. Minor deviations were observed in the underestimation of $d_{10}$ and in the slightly larger confidence intervals of $d_{43}$ at higher concentrations. These trends are consistent with the physical interpretation of the CFD-PBM model. In fact, increased particle number densities at higher concentrations naturally lead to higher collisions, broadening the PSD. A key outcome of this work is the demonstration that the proposed methodology enables, for the first time, the use of fully three-dimensional CFD-PBM simulations equipped with a rigorously calibrated kinetic model without resorting to dimensionality simplifications or ad hoc assumptions. This represents a substantial advancement over conventional parameter estimation workflows, which are typically limited by computational cost or require reduced-order approximations. By providing a modular, scalable, and inverse-design framework, this study opens new opportunities for the deployment of high-fidelity digital twin like models in precipitation, crystallization, and reactive mixing processes, where experimental characterization is limited, expensive, or intrinsically challenging. Overall, the results confirm that the hybrid CFD-PBM-DL framework constitutes a powerful and generalizable tool for identifying kinetic parameters in multivariate, nonlinear reactive systems, marking a significant step towards next-generation model-based process design and control. While this work focuses exclusively on a deep learning neural network as the mirror model, the primary objective was to establish a proof-of-concept for the inverse-design framework applied for the first time to 3D CFD-PBM simulations. Alternative machine learning algorithms, such as tree-based methods (Random Forest, XGBoost), could be explored in future studies to evaluate performance, robustness, and applicability to different datasets, but this was beyond the scope of the present study.

# GitHub

## Population Balance Model
Inverse design

## CRediT authorship contribution statement

**Antonello Raponi:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Formal analysis, Data curation, Conceptualization.
**Daniele Marchisio:** Writing – review & editing, Resources, Project administration, Funding acquisition.

# Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

# Acknowledgments

This project has received funding from ICSC – Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing (funded by the European Union – NextGenerationEU) and from PRIN (Progetti di ricerca di rilevante interesse nazionale – Bando 2022 – Project No. 2022JJRH8H: Non-equilibrium self-assembly of structured fluids: a multi-scale engineering problem). Computational resources were provided by HPC@POLITO, a project of Academic Computing within the Department of Control and Computer Engineering at the Politecnico di Torino (http://www.hpc.polito.it). This output reflects only the author's view. The European Health and Digital Executive Agency (HaDEA) and the European Commission cannot be held responsible for any use that may be made of the information contained therein.

# Data availability

Data will be made available on request.

# References

[1] T.E. Graedel, R. Barr, C. Chandler, T. Chase, J. Choi, L. Christoffersen, E. Friedlander, C. Henly, C. Jun, N.T. Nassar, D. Schechner, S. Warren, M. yu Yang, C. Zhu, Methodology of metal criticality determination, Environ. Sci. Technol. 46 (2012) 1063-1070.

[2] T.E. Graedel, E.M. Harper, N.T. Nassar, P. Nuss, B.K. Reck, Criticality of metals and metalloids, Proc. Natl. Acad. Sci. 112 (2015) 4257-4262.

[3] H. Béarat, M.J. McKelvy, A.V. Chizmeshya, R. Sharma, R.W. Carpenter, Magnesium hydroxide dehydroxylation/carbonation reaction processes: Implications for carbon dioxide mineral sequestration, J. Am. Ceram. Soc. 85 (2002) 742-748.

[4] X. Chen, J. Yu, S. Guo, Structure and properties of polypropylene composites filled with magnesium hydroxide, J. Appl. Polym. Sci. 102 (2006) 4943-4951.

[5] H. Gui, X. Zhang, W. Dong, Q. Wang, J. Gao, Z. Song, J. Lai, Y. Liu, F. Huang, J. Qiao, Flame retardant synergism of rubber and Mg(OH)2 in EVA composites, Polymer 48 (2007) 2537-2541.

[6] S. Zhang, F. Cheng, Z. Tao, F. Gao, J. Chen, Removal of nickel ions from wastewater by Mg(OH)2/MgO nanostructures embedded in Al2O3 membranes, J. Alloys Compd. 426 (2006) 281-285.

[7] H. Cao, H. Zheng, J. Yin, Y. Lu, S. Wu, X. Wu, B. Li, Mg(OH)2 complex nanostructures with superhydrophobicity and flame retardant effects, J. Phys. Chem. C 114 (2010) 17362-17368.

[8] C.M. Tai, R.K. Li, Studies on the impact fracture behaviour of flame retardant polymeric material, Mater. Des. 22 (2001) 15-19.

[9] S. Kakaraniya, C. Kari, R. Verma, A. Mehra, Gas absorption in slurries of fine particles: So2 - Mg(OH) 2 - MgSO3 system, Ind. Eng. Chem. Res. 46 (2007) 1904-1913.

[10] A. Sierra-Fernandez, L.S. Gomez-Villalba, O. Milosevic, R. Fort, M.E. Rabanal, Synthesis and morpho-structural characterization of nanostructured magnesium hydroxide obtained by a hydrothermal method, Ceram. Int. 40 (8 PART A) (2014) 12285-12292.

[11] X. Song, K. Tong, S. Sun, Z. Sun, J. Yu, Preparation and crystallization kinetics of micron-sized Mg(OH)2 in a mixed suspension mixed product removal crystallizer, Front. Chem. Sci. Eng. 7 (2013) 130-138.

[12] A. Cipollina, M. Bevacqua, P. Dolcimasco, L. Tamburini, A. Brucato, H. Glade, L. Buether, G. Micale, Reactive crystallisation process for magnesium recovery from concentrated brines, Desalination Water Treat. 55 (2015) 2377-2388.

[13] M. Shiea, A. Querio, A. Buffo, G. Boccardo, D. Marchisio, CFD-PBE modelling of continuous Ni-Mn-Co hydroxide co-precipitation for Li-ion batteries, Chem. Eng. Res. Des. 177 (2022) 461-472.

[14] T. Schikarski, M. Avila, H. Trzenschiok, A. Güldenpfennig, W. Peukert, Quantitative modeling of precipitation processes, Chem. Eng. J. 444 (2022) 136195.

[15] S. Romano, S. Trespi, R. Achermann, G. Battaglia, A. Raponi, D. Marchisio, M. Mazzotti, G. Micale, A. Cipollina, The role of operating conditions in the precipitation of magnesium hydroxide hexagonal platelets using NaOH solutions, Cryst. Growth Des. 23 (2023) 6491-6505.
[16] L. Lu, Q. Hua, J. Tang, Y. Liu, L. Liu, B. Wang, Reactive crystallization kinetics of magnesium hydroxide in the Mg(NO 3 ) 2 -NaOH system, Cryst. Res. Technol. 53 (2018) 1700130.

[17] J. Kennedy, R. Eberhart, Particle swarm optimization, in: Proceedings of ICNN'95-International Conference on Neural Networks, vol. 4, IEEE, 1995, pp. 1942–1948.

[18] F. Duan, F. Song, S. Chen, M. Khayatnezhad, N. Ghadimi, Model parameters identification of the PEMFCs using an improved design of crow search algorithm, Int. J. Hydrog. Energy 47 (2022) 33839–33849.

[19] A.G. Gad, Correction to: Particle swarm optimization algorithm and its ap-
plications: A systematic review, Arch. Comput. Methods Eng. 30 (2023)
3471–3471.

[20] Y. Chen, Z. Ye, B. Gao, Y. Wu, X. Yan, X. Liao, A robust adaptive hierarchical
learning crow search algorithm for feature selection, Electronics 12 (2023) 3123.

[21] V. Kumar, S.M. Yadav, A state-of-the-art review of heuristic and metaheuristic optimization techniques for the management of water resources, Water Supply 22 (2022) 3702–3728.

[22] B.A. Attea, A.D. Abbood, A.A. Hasan, C. Pizzuti, M. Al-Ani, S. Özdemir, R.D. Al-Dabbagh, A review of heuristics and metaheuristics for community detection in complex networks: Current usage, emerging development and future directions, Swarm Evol. Comput. 63 (2021) 100885.

[23] M. Prieß, J. Piwonski, S. Koziel, A. Oschlies, T. Slawig, Accelerated parameter identification in a 3D marine biogeochemical model using surrogate-based optimization, Ocean. Model. 68 (2013) 22–36.

[24] A. Raponi, D. Marchisio, Deep learning for kinetics parameters identification: A novel approach for multi-variate optimization, Chem. Eng. J. 489 (2024).

[25] I. Pan, L.R. Mason, O.K. Matar, Data-centric engineering: integrating simulation, machine learning and statistics. Challenges and opportunities, Chem. Eng. Sci. 249 (2022) 117271.

[26] A.I.J. Forrester, A. Sóbester, A.J. Keane, Engineering Design via Surrogate Modelling, Wiley, 2008.

[27] K. Jiang, Z. Liang, H. Jiang, Y. Luan, X. Su, T. Zheng, M. Liu, Y. Feng, W. Li, Y. Chen, G. Lu, Z. Liu, Systemic comparison of machine learning models in the optimization of flow field design for proton exchange membrane fuel cells, Energy 335 (2025) 138029.

[28] T. Chen, R. Li, X. Hu, B. Zhang, Y. Liu, L.L. Wang, N. Gao, Machine learning as CFD surrogate models for rapid prediction of building-related physical fields: A review of methods and state-of-the-art, Build. Environ. 285 (2025) 113667.

[29] X. Li, W. Sun, C. Qin, Y. Yan, L. Zhang, J. Tu, Evaluation of supervised machine learning regression models for CFD-based surrogate modelling in indoor airflow field reconstruction, Build. Environ. 267 (2025) 112173.

[30] L.-T. Zhu, X.-Z. Chen, B. Ouyang, W.-C. Yan, H. Lei, Z. Chen, Z.-H. Luo, Review of machine learning for hydrodynamics, transport, and reactions in multiphase flows and reactors, Ind. Eng. Chem. Res. 61 (2022) 9901–9949.

[31] D.L. Marchisio, J.T. Pikturna, R.O. Fox, R.D. Vigil, A.A. Barresi, Quadrature method of moments for population-balance equations, AIChE J. 49 (2003) 1266–1276.

[32] D.L. Marchisio, R.D. Vigil, R.O. Fox, Quadrature method of moments for aggregation-breakage processes, J. Colloid Interface Sci. 258 (2003) 322–334.

[33] A. Raponi, S. Romano, G. Battaglia, A. Buffo, M. Vanni, A. Cipollina, D. Marchisio, Computational modeling of magnesium hydroxide precipitation and kinetics parameters identification, Cryst. Growth Des. 23 (7) (2023) 4748–4759.

[34] A. Raponi, R. Achermann, S. Romano, S. Trespi, M. Mazzotti, A. Cipollina, A. Buffo, M. Vanni, D. Marchisio, Population balance modelling of magnesium hydroxide precipitation: Full validation on different reactor configurations, Chem. Eng. J. (2023) 146540.

[35] A. Raponi, D. Fida, F. Vicari, A. Cipollina, D. Marchisio, Computational fluid dynamics and population balance model enhances the smart manufacturing and performance optimization of an innovative precipitation reactor, Processes 13 (2025).

[36] D.L. Marchisio, R.O. Fox, A.A. Barresi, G. Baldi, On the comparison between presumed and full PDF methods for turbulent precipitation, 40, American Chemical Society, 2001, pp. 5132–5139,

[37] R. David, P. Marchal, J.P. Klein, J. Villermaux, Crystallization and precipitation engineering-III. A discrete formulation of the agglomeration rate of crystals in a crystallization process, Chem. Eng. Sci. 46 (1991) 205–213.

[38] D.L. Marchisio, A.A. Barresi, CFD simulation of mixing and reaction: The relevance of the micro- mixing model, Chem. Eng. Sci. 58 (2003) 3579–3587.

[39] D.L. Marchisio, A.A. Barresi, M. Garbero, Nucleation, growth, and agglomeration in barium sulfate turbulent precipitation, AIChE J. 48 (2002) 2039–2050.

[40] A. Mersmann, Crystallization Technology Handbook, Marcel Dekker, 2001, p. 832.