Application of Artificial Intelligence in Mathematical Modeling and Numerical Investigation of Transport Processes in Electromembrane Systems

1. Introduction

The modern era is characterized by the rapid development of artificial intelligence (AI) technologies, which are actively penetrating many areas of human activity, including membrane electrochemistry [1,2,3,4,5,6,7,8,9,10]. For example, Abuwatfa, AlSawaftah, et al. [1] present a systematic review of the application of artificial neural networks to predict membrane fouling in pressure-driven processes (microfiltration, ultrafiltration, nanofiltration, and reverse osmosis). The work demonstrates that artificial neural networks (ANNs), thanks to their nonlinear mapping capabilities, can effectively capture complex relationships between process parameters and fouling development, outperforming traditional empirical and mechanistic models. The article discusses the advantages and limitations of ANN approaches, identifying the potential of the method for a deeper understanding of fouling mechanisms.

In their work [2], Jawad, Hawari, and Zaidi developed an innovative hybrid approach combining artificial neural networks with response surface methodology (RSM) to model the forward osmosis process and predict permeate flux. This synergistic method not only captures nonlinear relationships between process parameters (solution concentration, temperature, membrane characteristics) but also performs sensitivity analysis to identify the most critical variables. The combined ANN-RSM approach provides both high predictive accuracy and a deep understanding of the mechanisms, enabling process optimization by identifying optimal parameter combinations.

Di Martino, Abraham, Pistikopoulos, and colleagues [3] applied rectified linear unit (ReLU) artificial neural networks to comprehensively optimize the performance of industrial reverse osmosis plants based on a superstructure approach. ANNs are used as surrogate models to quickly calculate the performance of various plant configurations, enabling engineers to perform global optimization while minimizing energy and capital costs. The model achieved exceptionally high accuracy in predicting plant current, salt drift, and other critical parameters, demonstrating the practical value of the approach for desalination plant design and operation.

In [4], Galinha and Crespo present a fundamental review of the evolution of machine learning in membrane process modeling over 25 years, tracing the transition from traditional “black box” approaches to modern interpretable approaches. The authors discuss various machine learning methods (ANN, PLS, PCA, chemometrics) applied to membrane performance modeling, fouling assessment, real-time process monitoring, and operational optimization. This article highlights a paradigm shift in the perception of machine learning—from an underrepresented scientific tool to a powerful method that, when properly applied, not only predicts but also improves our understanding of complex membrane processes.

The authors of [5] review current methods for predicting fouling in membrane bioreactors (MBRs), where membrane fouling remains a major barrier to the widespread use of this highly effective wastewater treatment technology. The paper critically analyzes the use of artificial neural networks and hybrid approaches for fouling prediction, discussing key challenges—including data quality, the risk of overfitting, and model interpretability issues. The article provides practical recommendations for the successful application of ANNs, emphasizing the need for combining multiple methods, long-term data collection, and continuous model improvement.

In [6], the authors present AI algorithms and a corresponding computing platform designed to predict membrane fouling processes in anaerobic membrane bioreactors (AnMBRs) used for wastewater treatment. The input factors included operating parameters, biomass properties, and membrane characteristics. The use of ANN and random forest methods made it possible to confirm the influence of individual factors and their interactions on the membrane fouling process, identifying important correlations necessary for optimizing their operation. The work [7] is devoted to the development of effective forecasting tools facilitating the implementation of organic solvent nanofiltration technology. An extensive dataset was compiled, including more than 18 characteristics, based on which AI models (neural networks, support vector machines, and random forest) were built to estimate the retention and transmission of molecules. Principal component analysis (PCA) showed the similarity of key factors determining both of these characteristics. The study [8] describes a comprehensive mathematical model for improving the performance of reverse osmosis systems in treating brackish water with moderate to high salt content. This model is based on the RSM and ANN approach, developed specifically for improving the performance of RO systems. The use of ANN showed significant superiority over the RSM model in terms of prediction accuracy. A two-stage optimization procedure made it possible to improve the system performance by achieving an optimal balance between the performance index and the energy consumption. The paper [9] discusses the use of a three-layer residual artificial neural network (R-TNN) for predicting the performance of thin-film nanocomposite reverse osmosis membranes (TFN-ROM). The resulting model demonstrated the ability to accurately account for the complex relationship between membrane permeability and retention capacity, which allows for the creation of highly efficient membrane structures. Finally, work [10] presents a deep learning protocol for predicting the electrical conductivity of hydroxide ions (OH-) in anion-exchange membranes (AEMs) based on poly (2,6-dimethylphenylene oxide) with functionalized ionic groups. The developed model achieves an accuracy of approximately 99,7%, which opens new prospects for the design of membranes with predetermined properties. This study emphasizes the effectiveness of an approach that consists of replacing standard physical variables with complex quantities constructed from these same quantities in specific combinations determined by the nature of the phenomena under study. This approach has several important advantages. First, a significant reduction in the number of independent variables facilitates the formulation and solution of the problem. Second, the use of complex variables that simultaneously integrate multiple influences of factors allows for a more detailed understanding of the internal structure and interrelations characteristic of the processes under study. As a result, the general form of the equations and graphical representation becomes more transparent and informative.

Another remarkable feature of complex variables is the fact that each fixed value of such a variable can correspond to an infinite number of different combinations of the basic quantities. Thus, each point in the space of generalized variables actually reflects an entire family of equivalent states united by common qualitative features. This fact allows for the simultaneous analysis of not just isolated particular cases, but entire classes of similar processes, increasing the level of generality of the conclusions obtained. The use of generalized dimensionless variables forms the basis of a specialized scientific field called similarity theory and dimensional analysis [11]. Although this approach has long been successfully used in physics, its adaptation to chemical disciplines remains an underdeveloped area of science. The goal of this work is to demonstrate the potential for applying similarity theory directly to membrane electrochemistry, revealing a fundamentally new perspective on traditional problems and opening broad horizons for a deeper understanding of the mechanisms of ongoing processes. As is well known, water is a vital resource for life, and the primary solution to the problem of freshwater shortage is access to effective desalination technologies, which include, in particular, electromembrane technologies. The efficiency of electromembrane systems for water purification, such as electrodialysis, is highly dependent on the hydrodynamics of the process, as well as the properties of the membrane surface. To reduce the costs of EMS experiments to intensify desalination, it is necessary to develop mathematical models of EMS water purification. Electroconvection is the primary mechanism of over-limit transport in membrane systems, and therefore, its intensification should be considered: an earlier onset to reduce the “plateau,” an increase in the size of electroconvective vortices, etc. When analyzing electroconvection processes, one must deal with a large number of diverse variables, such as the initial concentration of the substance, the change in electrical potential, the solution feed rate, the channel width, and other parameters, each of which acts as a separate factor. However, the influence of the specified parameters is not manifested individually, but in a complex, forming certain dimensionless quantities (criterion numbers, trivial similarity criteria such as Reynolds, Peclet numbers, etc., and the number of parameters decreases, according to the theorem of similarity theory. Criterion numbers characterize internal connections, for example, the Reynolds number characterizes the ratio of inertial forces to friction forces, and the Peclet number the ratio of convective transfer to diffusion, etc. and, depending on their value, determine the qualitative and quantitative picture of the solution flow and the transfer process. Moreover, the relationship between the physical parameters of the electroconvection process can be represented as a relationship between the similarity criteria, including such characteristics important for intensification as, for example, the dimensionless diameter of the electroconvective vortex (ECV), the length of the plateau of the current-voltage characteristic, ECV_CEM - the beginning of electroconvection near the CEM, ECV_AEM – the beginning of electroconvection near AEM.

1)

in the potentiostatic case we have (using the example of the dimensionless diameter of an electroconvective vortex D_ecv):

$D_{ecv}=f(\mathrm{Re},Pe,\epsilon ,K_{el},d_{\phi })$, where $d_{\phi }$ – dimensionless potential jump;

2) in the potentiodynamic case we have:

$D_{ecv}=f(\mathrm{Re},Pe,\epsilon ,K_{el},d_1)$, where $d_1$ – dimensionless potential sweep rate. Besides, Re – Reynolds number, Pe — Peclet number, $\epsilon$ – this is the ratio of the square of the thickness of the equilibrium space charge region to the square of the distance between the membranes [12,13], and $K_{el}$ – this is the ratio of the electrical force to the inertial force. Typically, to construct functions like $f$ either an exact solution to the boundary value problem of a mathematical model or experimental studies are used. Currently, finding an analytical solution to the electroconvection problem is difficult due to insurmountable mathematical difficulties, and an experimental determination is extremely labor-intensive due to the large number of criterion numbers. In this regard, this paper proposes using neural networks to determine functional dependencies. To implement neural networks, a sample was compiled based on numerical experiments with the developed mathematical model. The entire sample was divided into two sets: a set used for training (setting weights and biases) and a test set, which was used to determine the correctness of the neural network. A multilayer feedforward network was chosen as the network architecture; the backpropagation algorithm with different optimizers was used to train the neural network. The application of the found functions can contribute to the development of a methodology for calculating the optimal geometric and technological parameters of the desalination process in electromembrane units.

2. Methods

2.1. Mathematical Model

The flow of liquid in a rectangular channel of length L and width H, equipped with cation-exchange and anion-exchange membranes, is studied; a schematic representation is shown in Figure 1. Modeling the dynamics of liquid flow in a flow channel, taking into account the effects of electroconvection, is carried out by means of a joint solution of the system of Nernst-Planck-Poisson equations (formulas 1-4) and the Navier-Stokes equations (formulas 5-6).

$${\overrightarrow{j}}_i=-{\displaystyle \frac{F}{RT}}{z}_i{D}_i{C}_i\overrightarrow{E}-D_i\nabla C_i+C_i\overrightarrow{V}, i=1,2$$

(1)

$${\displaystyle \frac{\partial C_i}{\partial t}}=-div{\overrightarrow{j}}_i, i=1,2$$

(2)

$${\epsilon }_r\Delta \phi =-F(z_1{C}_1+z_2{C}_2)$$

(3)

$$\overrightarrow{I}=F(z_1{\overrightarrow{j}}_1+z_2{\overrightarrow{j}}_2)$$

(4)

$${\displaystyle \frac{\partial \overrightarrow{V}}{\partial t}}+(\overrightarrow{V}\nabla )\overrightarrow{V}=-{\displaystyle \frac{1}{{\rho }_0}}\nabla P+\nu \Delta \overrightarrow{V}+{\displaystyle \frac{1}{{\rho }_0}}\overrightarrow{f}$$

(5)

$$div(\overrightarrow{V})=0$$

(6)

where

$$\overrightarrow{f}=\rho \overrightarrow{E}=-{\epsilon }_r\Delta \phi \overrightarrow{E}={\epsilon }_r\Delta \phi \nabla \phi ={\epsilon }_r\overrightarrow{E}div\overrightarrow{E}$$

.

Here $\overrightarrow{V}$– flow velocity of the solution; Index 1 refers to cations, index 2 refers to anions; ${\overrightarrow{j}}_1, {\overrightarrow{j}}_2$ – fluxes; $C_1, C_2$ – concentrations in the solution;$z_1, z_2$ – charge numbers; $\overrightarrow{I}$ – current density; $D_1,D_2$ – diffusion coefficients; $\phi$ – electric potential field; $\overrightarrow{E}=-\nabla \phi$ – electric field strength; ${\epsilon }_r$ – dielectric permittivity of the electrolyte; $F$ – Faraday’s constant; $R$ – universal gas constant; $T$ – absolute temperature; $t$– time; ${\rho }_0$ – density; $\nu$ – kinematic viscosity; $P$ – pressure; $\rho =F(z_1{C}_1+z_2{C}_2)$ – distribution of charge densities.

The potentiodynamic regime is being investigated. The boundary conditions are set within the following constraints:

1) Boundary conditions on the membrane surfaces

${С}_2\left(t,0,y\right)={С}_{2m}$, $\left({\displaystyle \frac{\partial C_1}{\partial x}}-{\displaystyle \frac{F}{RT}}{С}_{1}E\right)(t,0,y)=0$, $\phi \left(t,0,y\right)=0$

${С}_1\left(t,H,y\right)=C_{1m}(t)$, $\left({\displaystyle \frac{\partial C_2}{\partial x}}+{\displaystyle \frac{F}{RT}}{С}_{2}E\right)(t,H,y)=0$, $\phi (t,H,y)=\Delta \phi (t)$

${С}_1\left(0,x,y\right)={С}_0$, ${С}_2\left(0,x,y\right)={С}_0$, $\phi \left(0,x,y\right)=0$,

where $\Delta \phi (t)$ denotes a potential jump similar to that indicated earlier, and is traditionally performed $\Delta \phi (t)=-d_1\cdot t$. The rate of increase in the potential jump is determined by the parameter d₁, whose unit is volts per second [V/s].

The no-slip condition is assumed for the fluid velocity on the channel walls.

2) Boundary conditions at the channel inlet: At the inlet to the desalination channel, specified ion concentrations are established that satisfy the condition of electroneutrality, that is:

$$C_i(t,x,0)=C_{i,0}i=1,2$$

$$\phi (t,x,0)=\Delta \phi \cdot x/H$$

(7)

A constant concentration distribution, coupled with Ohm’s law, leads to the expression presented in formula (7). Consequently, this formula is consistent with both Ohm’s law and the condition of constant concentrations at the channel entrance.

Distribution of solution flow rate $\overrightarrow{V}$ at the input we accept the corresponding Poiseuille law, according to which the velocity profile has a parabolic form, but $V_0$ denotes the average flow rate of the solution.

$$V_x=0, V_y=6V_0{\displaystyle \frac{x}{H}}(1-{\displaystyle \frac{x}{H}}).$$

(8)

3) Boundary conditions at the outlet of the desalination channel (DC):

a) At the boundary of the output region $y=L, x\in [0,H], t\ge 0$ for component concentrations, a condition ensuring conservation of ion flux is imposed, meaning that removal of salt ions from the DC area happens solely due to convective transport:

$-\overrightarrow{n}\cdot {\overrightarrow{j}}_i(x,L,t)=-V_y(x,L,t)C_i(x,L,t), i=1,2$

b) For the potential jump the following condition is set:

$-\overrightarrow{n}\cdot \nabla \phi =0$,

i.e.

${\displaystyle \frac{\partial \phi (x,L,t)}{\partial y}}=0$

4) The initial conditions are set as follows: For concentrations and electric potential:

$C_i(0,x,y)=C_{i,0}i=1,2$.

$\phi (0,x,y)=\Delta \phi \cdot x/H$

We assume that the flow velocity of the liquid is distributed according to Poiseuille’s law (8).

2.2. Similarity Method for Transport Processes in a Flat Channel

This work continues and develops the ideas presented in studies [14,15,16], where a theory of process similarity was developed for a two-dimensional desalination channel.

1. Conversion to dimensionless form

We select characteristic parameters as the basic scaling quantities $C_0,V_0,d_{\phi },H$ and we will make the transition to a dimensionless representation, considering the following expressions to be valid:

$x^{(u)}={\displaystyle \frac{x}{H}};$ $y^{(u)}={\displaystyle \frac{y}{H}};$ $L^{(u)}={\displaystyle \frac{L}{H}};$ $t^{(u)}={\displaystyle \frac{t{V}_0}{H}};$ ${\overrightarrow{V}}^{(u)}={\displaystyle \frac{\overrightarrow{V}}{V_0}};$ $P^{(u)}={\displaystyle \frac{P}{{\rho }_0{V}_0}};$

${\overrightarrow{E}}^{(u)}={\displaystyle \frac{HF}{R{T}_0}}\overrightarrow{E};$$D_i^{(u)}={\displaystyle \frac{D_i}{D_0}},$$i_{av}{}^{\left(u\right)}={\displaystyle \frac{i_{av}H}{D{C}_{0}F}};$ ${\eta }^{\left(u\right)}={\displaystyle \frac{\eta }{D{C}_{0}F}};$ $f_0={\displaystyle \frac{{\rho }_0{V}_0^2}{H}};$

$d_{\varphi }^{(u)}={\displaystyle \frac{F}{R{T}_0}}{d}_{\varphi };$${\phi }^{(u)}={\displaystyle \frac{F}{R{T}_0}}\phi ;$ ${\epsilon }^{(u)}={\displaystyle \frac{R{T}_0\epsilon }{F^2{H}^2{C}_0}};$ ${С}_i^{(u)}={\displaystyle \frac{C_i}{C_0}};$

${\overrightarrow{j}}_i^{(u)}={\displaystyle \frac{H}{D_0{C}_0}}{\overrightarrow{j}}_i;$${\overrightarrow{I}}^{(u)}={\displaystyle \frac{H}{D_0{C}_{0}F}}\overrightarrow{I};$

where $D_0=D_1$ – average diffusion coefficient, and the superscript (u) denotes that the quantity under consideration is dimensionless.

2. Physical meaning of dimensionless parameters in equations and boundary conditions

In paper [15], the boundary value problem associated with model [14] was reformulated into a dimensionless form by utilizing the characteristic scaling factors mentioned earlier. This transformation led to the introduction of several key dimensionless parameters:

$L^{(u)}={\displaystyle \frac{L}{H}}$ – representing the dimensionless channel length, $Pe={\displaystyle \frac{V_{0}H}{D_0}}$ – known as the Peclet number, indicating the ratio of convective transport intensity to molecular diffusion, $\mathrm{Re}={\displaystyle \frac{V_{0}H}{\nu }}$ – referred to as the Reynolds number, illustrating the equilibrium between inertial $F_{in}={\rho }_0{H}^2{V}_0^2$ and viscous friction forces $F_{tr}=\nu {\rho }_0{V}_{0}H$, along with additional dimensionless numbers ${\epsilon }^{(u)}={\displaystyle \frac{RT\epsilon }{H^2{C}_0{F}^2}}$ and $K_{el}={\displaystyle \frac{RT{C}_0}{{\rho }_0{V}_0^2}}$.

3. Physical interpretation of the quantity ${\epsilon }^{(u)}$This parameter, first introduced in [12,13], has received the following physical interpretation. The formula can be represented as follows: ${\epsilon }^{(u)}={\displaystyle \frac{RT\epsilon }{H^2{C}_0{F}^2}}=2{\left[{\displaystyle \frac{l_d}{H}}\right]}^2$ , where $l_d=\sqrt{{\displaystyle \frac{RT\epsilon }{2C_0{F}^2}}}$– Debye length, which characterizes the size of the equilibrium space charge region.

Thus, the parameter ${\epsilon }^{(u)}$ is the ratio of the square of twice the Debye length (the thickness of the space charge region) to the square of the intermembrane distance, or equivalently, the ratio of the Debye length itself to half the intermembrane distance. Physically, this means that the parameter ${\epsilon }^{(u)}$ determines the ratio of the sizes of the charged and uncharged zones of the solution within the intermembrane gap.

4. Physical interpretation of the quantity $K_{el}={\displaystyle \frac{RT{C}_0}{{\rho }_0{V}_0^2}}$When the fluid flows with an average velocity $V_0$, the quantity $P_0={\rho }_0{V}_0^2$ serves as a characteristic pressure, often referred to as the velocity head or dynamic pressure. It is frequently utilized in both hydrodynamics and aerodynamics as a scaling factor for hydrodynamic pressures and appears in various aerodynamic coefficients.

The electric pressure acting on a unit cross-sectional area is written as:

$$P_{el}=2RT{C}_0$$

In addition, the value $RT{C}_0$, according to Van’t Hoff’s law corresponds to the osmotic pressure of a solution with a concentration of $C_0$.

Based on this, the value of the ratio $K_{el}={\displaystyle \frac{P_{el}}{P_0}}$ can be interpreted as the ratio of the characteristic pressure caused by the action of electrical force to the characteristic hydrodynamic pressure arising as a result of the movement of the solution.

5. The electroconvection criterion number $K_{ek}$ is characterized as the coefficient of the dimensionless electric force driving the generation of electroconvective flows. Its expression in terms of physical parameters is as follows:

$$K{}_{ek}={\epsilon }^{(u)}{K}_{el}.$$

Then, expressing both components through dimensional parameters, we obtain:

$$K{}_{ek}=2{\left[{\displaystyle \frac{l_d}{H}}\right]}^2{\displaystyle \frac{RT{C}_0}{{\rho }_0{V}_0^2}}=2{\displaystyle \frac{RT{C}_0{l}_d^2}{{\rho }_0{V}_0^2{H}^2}}=2{\displaystyle \frac{F_{el}}{F_{in}}},$$

where $F_{el}$ — electric force acting on a space charge region, and $F_{in}$ — force of inertia.

This expression shows that the criterion number of electroconvection $K_{ek}$ is the ratio of electrical force to inertial force.

It is important to note that the inertial force does not depend on the initial concentration C₀, since it is determined by the kinetic properties of the liquid. On the other hand, the electric force $F_{el}$, being the product of the electric pressure and the space charge volume, also turns out to be independent of the initial concentration. Increasing the concentration increases the electric pressure but simultaneously narrows the space charge volume, which mutually compensates for the effect of the initial concentration.

$$F{}_{el}=RT{C}_0{l}_d^2=RT{C}_0{\displaystyle \frac{RT\epsilon }{2C_0{F}^2}}={\displaystyle \frac{{(RT)}^2\epsilon }{2F^2}}$$

Therefore, the criterion number of electroconvection $K_{ek}$ does not actually depend explicitly on the initial concentration, which confirms the assumption about the weak connection between electroconvection and concentration.

The definition of the number $K_{ek}$ implies the presence of a certain critical value ${\overline{K}}_{ek}$, upon exceeding which the widespread formation of electroconvective structures in the channel begins.

An assessment of the dimensionless parameters reveals that for typical values of dimensional quantities relevant to electrodialysis, the Peclet and Reynolds numbers fall within the orders of magnitude of approximately ${10}^2÷{10}^4$ and $1÷100$ respectively. Additionally, the parameter $\epsilon$ spans a range from ${10}^{-17}$ to ${10}^{-7}$, making it possible to treat it as a small parameter.

The number $K_{el}$ is of the order of $10÷{10}^3$, and the number $K_{ek}$ is of the order of ${10}^{-14}÷{10}^{-4}$.

6. In addition to the dimensionless parameters considered earlier, additional dimensionless quantities are introduced, appearing in the boundary conditions [15]:

$d_{\phi }^{(u)}={\displaystyle \frac{F}{R{T}_0}}{d}_{\phi }$ – total potential value;

$d_0^{(u)}={\displaystyle \frac{F}{R{T}_0}}{d}_0$ – initial value of potential;

$d_1^{(u)}={\displaystyle \frac{F}{R{T}_0}}{d}_1\cdot t$ – growth rate of potential;

${С}_{km}^{(u)}={\displaystyle \frac{C_{km}}{C_0}}$ – limit value of concentration on a CEM;

${С}_{am}^{(u)}={\displaystyle \frac{C_{am}}{C_0}}$ – limit value of concentration on the AEM.

These parameters reflect the physical conditions present at the boundaries of the channel and guarantee proper formulation of the boundary value problem. It is postulated that the initial concentration within the channel matches the concentration at the channel inlet. For the sake of simplification in subsequent computational experiments, a unified initial concentration equivalent to the entry concentration has been implemented. $C_{am}=C_{km}=C_0$.

7. Dimensionless equations and boundary conditions

Upon transitioning to dimensionless variables and parameters, the system of Nernst – Planck – Poisson and Navier – Stokes equations, originally stated in [14], takes the following form (for notational convenience, the dimensionless index “u” is dropped throughout):

$${\overrightarrow{J}}_i=z_i{D}_i{C}_i\overrightarrow{E}-D_i\nabla C_i+Pe{C}_i\overrightarrow{V},i=1,2$$

$$Pe{\displaystyle \frac{\partial C_i}{\partial t}}=-div{\stackrel{\rightharpoonup }{j}}_i,i=1,2$$

$$\epsilon \Delta \phi =-\left(z_1{C}_1+z_2{C}_2\right)$$

$$\overrightarrow{I}=z_1{\overrightarrow{j}}_1+z_2{\overrightarrow{j}}_2$$

$${\displaystyle \frac{\partial \overrightarrow{V}}{\partial t}}+(\overrightarrow{V}\nabla )\overrightarrow{V}=-\nabla P+{\displaystyle \frac{1}{\mathrm{Re}}}\Delta \overrightarrow{V}+\epsilon K_{el}\Delta \phi \nabla \phi$$

$$div(\overrightarrow{V})=0$$

The Nernst–Planck–Poisson and Navier–Stokes equations embody expressions of basic conservation principles. Consequently, the comprehensive model describing ion transport incorporating electroconvection, beyond merely defining particular equations, primarily hinges on the specification of appropriate boundary conditions. Detailed information regarding dimensional boundary conditions is available in [14]. These include:

$x=0$– the boundary between the surface of the AEM and the solution;

$x=1$ – the boundary between the surface of the CEM and the solution;

$y=0$ – inlet to the desalination channel;

$y=L$ – outlet from the desalination channel.

Let us now write the dimensionless boundary conditions considering the previously introduced dimensionless parameters.

1. Conditions on Ion Exchange Membrane Surfaces

a) Surface of Anion Exchange Membrane: ($x=0,y\in [0,L],t\ge 0$):

$${С}_2(t,0,y)={С}_{am}$$

Additionally, we presume that the AEM exhibits ideal selectivity and equipotentiality characteristics:

$$\left({\displaystyle \frac{\partial C_1}{\partial x}}+z_1{C}_1{\displaystyle \frac{\partial \phi }{\partial x}}\right)(t,0,y)=0.$$

$$\phi (t,0,y)=0$$

b) Surface of Cation Exchange Membrane: ($x=1,y\in [0,L],t\ge 0$):

$${С}_1(t,1,y)={С}_{km}$$

$$\left({\displaystyle \frac{\partial C_2}{\partial x}}+z_2{C}_2{\displaystyle \frac{\partial \phi }{\partial x}}\right)(t,1,y)=0.$$

$$\phi (t,1,y)=d_{\phi }(t),$$

where $d_{\phi }=d_0+d_1\cdot t$, $t\ge 0,y\in [0,L]$, $d_0,d_1$ – the initial value and the growth rate of the potential jump, respectively.

A no-slip condition is imposed on the membrane surface, meaning zero velocity of the liquid relative to the walls:

$$\begin{array}{l}{V}_x(t,0,y)=0,V_y(t,0,y)=0,\\ V_x(t,1,y)=0,V_y(t,1,y)=0\end{array}$$

2. At the inlet to the DC $x\in \left[\mathrm{0,1}\right],y=0,t\ge 0$At the entrance to the channel, the concentrations of sodium and chlorine ions are assumed to be uniformly distributed, and the condition of electrical neutrality is assumed to be met, that is, the equality of positive and negative charges.

${С}_i(t,x,0)=1,i=1,2$.

The condition is imposed on the potential:

${\displaystyle \frac{\partial \phi (t,x,0)}{\partial y}}=-{\displaystyle \frac{1}{z_1^2{D}_1+z_2^2{D}_2}}\left(z_1{D}_1{\displaystyle \frac{\partial C_1(t,x,0)}{\partial y}}\right.\left.+z_2{D}_2{\displaystyle \frac{\partial C_2(t,x,0)}{\partial y}}\right)$

ensuring the absence of currents.

The velocity profile at the channel entrance has a parabolic shape corresponding to the Poiseuille distribution.

$V_x\left(t,x,0\right)=0,$ $V_y\left(t,x,0\right)=6x\left(1-x\right)$

3. At the outlet of the DC $x\in \left[\mathrm{0,1}\right],y=L,t\ge 0$

For concentrations, we’ll apply the condition for ion flow at the outlet, assuming that the removal of salt ions from the DC results purely from the flow of the solution itself:

$z_i{C}_i(t,x,L){\displaystyle \frac{\partial \phi (t,x,L)}{\partial y}}+{\displaystyle \frac{\partial C_i(t,x,L)}{\partial y}}=0i=1,2$

${\displaystyle \frac{\partial \phi (t,x,L)}{\partial y}}=0$

Velocity conditions at the inlet and outlet are alike:

$$V_x(x,L,t)=0,V_y(x,\mathrm{L},t)=6x\left(1-x\right)$$

.

Moreover, at the corners of the channel’s outlet, an extra pressure condition is applied, permitting free departure of vortex formations from the desalination channel domain.

4. Initial Conditions

Initial conditions ($t=0$) conform to the given boundary conditions. As an illustration, the initial ion concentration inside the channel aligns with the concentration level at the channel inlet, and for the potential $\phi$ the initial condition can be a simple linear function independent of $y$:

$$\phi (0,x,y)=d_{\phi }x$$

Such conditions ensure a smooth start of the modeling process and contribute to rapid stabilization of the numerical solution.

According to the above reasoning, electroconvection depends not explicitly, but indirectly on the initial concentration $C_0$ (and boundary concentrations $C_{am}$, $C_{km}$).

It has been experimentally established (Figure 2) that increasing the boundary concentrations $C_{am}=C_{km}=C_0$ by a factor of 10 results in only a slight change in the size of electroconvective vortices, measured along the outer contour of the ion trajectory envelope. This is due to the fact that the main factors determining the intensity of electroconvection are the hydrodynamic characteristics of the flow and the electric field, and not the concentration of the salt solution itself.

The influence of boundary concentrations $C_{am}=C_{km}=C_0$ on the current-voltage characteristic is also insignificant (Figure 3).

The protrusion (local maximum) between the ohmic region and the plateau in Figure 3a,b is explained by the non-stationarity effects caused by the high sweep speed d₁=0.49 (see section 3.1).

3. Results

3.1. Training Set for Constructing a Neural Network

Important characteristics of the electroconvective transfer process include the dimensionless diameter of the ECV, the length of the CVC plateau, the onset of electroconvection near the AEM, and the onset of electroconvection near the CEM. These characteristics are determined through a computational experiment, followed by recording of the corresponding values. However, conducting such computational experiments is time-consuming and resource-intensive. For example, conducting similar experiments in this study took from several days to several weeks, depending on the input parameters. Therefore, for such problems, it is advisable to learn how to predict the required characteristics while avoiding direct calculations. Neural networks can be used for this purpose [17], after first collecting a representative training sample through direct calculations.

Since the initial concentration has little effect on the sizes of the ECV and CVC, when conducting computational experiments to collect information for the training sample, it was decided to vary the potential sweep rate and the solution pumping rate, and the corresponding criterion numbers will also change.

Thus, the task is to find the following functional dependencies:

$D_{ecv}=f(\mathrm{Re},Pe,K_{el},d_1)$ – dimensionless diameter of an electroconvective vortex;

$L_{pl}=p(\mathrm{Re},Pe,K_{el},d_1)$ – length of the plateau of the current voltage characteristic;

$EC{V}_{AEM}=h(\mathrm{Re},Pe,K_{el},d_1)$ – the beginning of electroconvection near the AEM;

$EC{V}_{CEM}=g(\mathrm{Re},Pe,K_{el},d_1)$ – the beginning of electroconvection near the CEM, here $d_1$ – dimensionless potential sweep rate.

Let’s first vary the potential sweep rate, fix the initial concentration in $C_0=0,01$ mol/m³ and the initial velocity in $V_0={10}^{-4}$ m/s. We will determine the ECV diameter, the CVC plateau length, and the onset of electroconvection near the AEM and CEM, first in dimensional form, and then nondimensionalize using the formulas given in Section 2.2 of this study. To determine the plateau length, it is necessary to plot the CVC characteristic. Figure 4 shows the current-voltage characteristic at a potential sweep rate of 0,00375 V/s. L_pl was determined from region №2 of the CVC characteristic, and the ECV_CEM and ECV_AEM parameters from regions 3 and 4, respectively.

The onset of electroconvection on the membranes was also confirmed by the graph of the solution flow lines, and the diameter of the ECV was also determined from it (Figure 5).

As a result of conducting a series of computational experiments with different potential sweep rates and calculating the required characteristics in dimensionless form, we obtained Table 1.

Let us plot graphs of the dependence of the found parameters on the dimensionless potential sweep rate (Figure 6). Figure 6a suggests that the dimensionless diameter of the ECV is virtually independent of the potential sweep rate and varies from 0,32 to 0,36. The length of the CVC plateau is virtually constant at low sweep rates, but then decreases significantly at high d₁ (Figure 6b). The parameters of the onset of electroconvection on the CEM and AEM behave stably at low potential sweep rates, but at high values, their character becomes abrupt (Figure 6c,d). Thus, it can be concluded that a quasi-stationary regime is realized at low sweep rates (0.01 to 0.07), beyond which the desalination process becomes non-stationary. One manifestation of these nonstationarity effects, as seen by comparing the CVC in Figure 3 and Figure 4, is the appearance of a small hump near the limiting current on the CVC.

Table 1 was then expanded to Table 2 by conducting a series of experiments with the initial velocity and potential sweep rate at a constant initial concentration of $C_0=0,01$ mol/m³. The data from this table served as the training set for the neural network being developed.

3.2. Building a Neural Network

The neural network being developed should solve the task of simultaneously predicting four parameters (multi-output regression); however, before feeding the training sample (Table 2) into the neural network, preprocessing of the data must be performed.

Data normalization plays an important role in the training process of neural networks. It helps significantly increase training efficiency and improve the quality of the resulting models. Normalization is necessary for:

1. Optimizing training speed

Neural networks use the gradient descent method to minimize the objective function. Gradients are calculated relative to changes in weights, and if the input data have different scales, this slows down the training process. The greater the difference in feature scales, the longer it will take for the model to converge.

2. Improving accuracy and robustness

Most algorithms depend on the distance between features, and differences in scale can distort the distance between elements, negatively impacting model quality.

3. Addressing multicollinearity

Some features can be correlated with each other, degrading model performance.

Given the above reasons, the training set was normalized. After this, a data split ratio was chosen: 80% for training, 20% for testing, and cross-validation was performed into five parts. Cross-validation is necessary due to the small data volume, the need for an accurate final estimate, and when comparing multiple neural network parameter configurations for objective comparison. When selecting neural network parameters such as the number of layers, neurons, learning rate, activation functions, and optimizers, it is useful to use hyperparameter selection using Grid Search or Randomized Search. Since searching for the best hyperparameters by randomly selecting them from a given distribution of possible values is more efficient than Grid Search for large hyperparameter spaces, Randomized Search was chosen. The standard MSE (Mean Square Error) loss function for regression problems was used, training lasted 50 epochs. As a result of hyperparameter optimization, the best performance was achieved using an architecture of 4-3-4, the LBFGS optimizer, and logistic activation function. The LBFGS optimizer (Limited-memory Broyden-Fletcher-Goldfarb-Shanno) is often used when high convergence speed and good accuracy are required, especially for small or medium-sized datasets. Its main advantage lies in its ability to converge quickly even on smaller samples and efficient operation with limited memory usage. Graphs illustrating loss functions for the testing and training sets are shown in Figure 7, while predicted versus actual values for each forecast variable along with their corresponding root mean square errors (RMSE) are presented in Figure 8.

The following metrics were used to evaluate the final quality of the developed neural network: RMSE (Root Mean Square Error) for individual measurements, Total RMSE (the overall metric for all parameters), MAPE (Mean Average Percentage Error) for relative error, and R² Score (the coefficient of determination)—a measure of the correspondence between the predicted and true relationships. All obtained metrics are presented in Table 3.

We will verify the performance of the neural network on new data not included in either the training or test set, and calculate the percentage error (Table 4).

4. Discussion

This study proposes an effective hybrid approach combining neural networks with numerical models for analyzing electromembrane desalination processes, the key mechanism of which is electroconvection.

A multilayer feedforward neural network was developed, tuned using a randomized hyperparameter search, which improved the efficiency and accuracy of predicting system characteristics. The model achieved the following performance metrics: RMSE=0.417, MAPE=3.74% and R²=0.961. This indicates a strong agreement between the calculated and empirically observed parameters, emphasizing the importance of the proposed approach. The resulting neural network model demonstrated high accuracy in reproducing key process indicators—such as the dimensionless electroconvective vortex diameter, the current-voltage characteristic plateau length, and the onset of electroconvection near the CEM and AEM—based on input parameters including the potential sweep rate, Péclet number, Reynolds number, and $K_{el}$. The application of the proposed method opens broad prospects for membrane design and the development of adaptive desalination process control strategies.

A key challenge in potentiodynamic modeling is selecting an appropriate potential sweep rate. While a lower rate enhances accuracy, it proportionally increases computational expense. Conversely, a higher rate reduces computation time at the cost of introducing errors due to non-stationary effects. This study identifies a dimensionless sweep rate of 0.07 as the optimal threshold, representing the maximum value for which a quasi-steady-state regime is maintained, thus balancing computational efficiency with model fidelity.

This work represents an important contribution to the development of interdisciplinary research at the intersection of membrane electrochemistry, artificial intelligence and computer modeling.