Multiphysics Modeling of the Electrochemical Response of Screen-Printed Electrodes for Sensing Applications

In this work, we present and validate a new 3-D COMSOL Multiphysics model capable of simulating the cyclic voltammetry (CV) and the electrochemical impedance spectroscopy (EIS) response of screen-printed devices. The proposed model considers the dominant electrical and electrochemical phenomena on the electrodes and in the electrolyte solution with redox reactions occurring at the metal/solution interface. The terminals of the simulated device are virtually connected to an equivalent circuit of a potentiostat in order to apply the potential and perform the CV readings with possibility to extend to other types of voltametric measurements. The model is calibrated on screen-printed devices through dedicated sets of experimental measurements that use CV and EIS with 10 mM [Fe(CN)6] $^{{3}-/{4}-}$ as redox mediator in phosphate-buffered saline (PBS) solution. The model parameters for each device are retrieved from the experimental measurements, e.g., double layer capacitance, red/ox diffusivity, electrolyte conductivity equilibrium, and so on. After the calibration, the model is capable of simulating CV responses using different scan rates and redox couple concentrations. The obtained simulated CV signals are significantly consistent with the experimental data, fitting well the experimental curves of both the devices. The simulated/experimental peaks differ less than 50 mV, while the anodic and cathodic current peaks vary only of few microamperes ( $\lt 11~\mu $ A), with values within the standard deviation of measurement experimental measurements. Finally, the proposed model will allow the simulation of more complex electrode morphologies and electrochemical conditions as well, as the model parameters are retrieved experimentally through specific measurements.


I. INTRODUCTION
I N RECENT years, electrochemical biosensors have attracted considerable interest for their highly sensitive and Giulio Rosati is with CSIC and BIST, Catalan Institute of Nanoscience and Nanotechnology (ICN2), 08193 Barcelona, Spain.
Arben Merkoçi is with CSIC and BIST, Catalan Institute of Nanoscience and Nanotechnology (ICN2), 08193 Barcelona, Spain, and also with the Catalan Institution for Research and Advanced Studies (ICREA), 08010 Barcelona, Spain.
Digital Object Identifier 10.1109/JSEN.2024.3410094 fast detection performance and for their lower price [1], [2], [3].Due to their characteristics, these devices have attracted attention, especially for their application in smart healthcare and clinical diagnostics [4], [5], [6], [7], sustainable food production, and agriculture 4.0 [8], [9], [10], [11].While having already proved to be powerful resources exploitable in a wide range of sectors, the design of biosensor devices is often based on empirical evaluations of their performance without many insights into their modeling before the fabrication.Nowadays, the mathematical modeling of physical systems through finite-element simulations is essential for engineering and nanotechnology.Mathematical modeling allows the investigation of different designs and parametrical configurations of the system, evaluating its predicted behavior and outcome by changing, e.g., geometries, surface morphologies, proprieties, and testing conditions.
The commercial finite-element software COMSOL Multiphysics 1 is suitable for multiphysics simulations of a wide 1 Registered trademark.range of nonstandard geometries and physical phenomena, e.g., hydrodynamic, and electrochemical processes.However, only a few works report sensor modeling, and they often lack validation through experimental measurements [12], [13].In [14], electrochemical simulations in COMSOL Multiphysics focus on the standard approximations of electrochemical theory implemented in the software and review the concept of finite-element modeling of the complex electrochemical phenomena.Modeling the electrical responses of biosensors, such as cyclic voltammetry (CV), is complex due to the need: 1) to integrate synergetic multiphysics models; 2) to properly calibrate parameters related to device materials; and 3) to appropriately define the testing condition, with specific equivalent circuits for signal generation and conditioning, e.g., potentiostats for CVs and frequency response analyzers for EISs.
As reported in the literature, finite-element simulations in COMSOL Multiphysics are exploited to design optimized electrodes to enhance their performance [15], [16], [17].However, the proposed models seldom consider the full electrode characterization and investigate the physical, chemical, and electrical phenomena on (typically a small portion of) its surface.For instance, in [18], a proposed biosensor was modeled through COMSOL Multiphysics to study the o-quinone enzymatic production and its electrochemical behavior, with the results subsequently validated experimentally.Another study investigated the influence of the working electrode (WE) area and the working/counter electrode (CE) distance on the sensor's electrical performance [15].However, these models lacked a 3-D structure, and complex phenomena contributing to the electrochemical response were not fully considered.With our previous work [19], the electrochemical impedance spectroscopy (EIS) response of screen-printed electrodes (SPEs) was simulated in 3-D by a Multiphysics modeling approach, also with self-assembled monolayers of 11-mercaptoundecanoic acid (MUA).The study was conducted in the frequency domain, while time-dependent modeling for the simulation of potentiometric techniques remains unexplored.
In this work, we report a novel comprehensive multiphysics model for the simulation of CV of electrochemical sensors using COMSOL Multiphysics.The model is designed to account for all the relevant chemical and electrical phenomena characterizing the CV response of an electrochemical sensor with a redox couple in an electrolyte solution.First, we calibrated the model with experimental measurements performed on two commercial screen-printed devices with different electrode surface proprieties.Then, we validated the simulated CV responses by comparing with new experimental datasets obtained under different measurement conditions, i.e., different concentrations and scan rates.Our results constitute an essential starting point for finite-element modeling of electrochemical sensors, which is fundamental for future analysis of the electrochemical and electrical processes of sensors in order to study their characteristics, identify eventual issues, and optimize their performance.The developed model will allow as well the simulation of more complex electrode morphologies and electrochemical conditions, considering that the model parameters are not set with theo- retical values but with the ones obtained by the experimental measurements.

II. EXPERIMENTAL DEVICES AND MEASUREMENT SETUP A. Devices
The electrochemical devices under test were commercial SPEs from Metrohm DropSens, Spain.The devices were fabricated on 34 × 10 × 0.5 mm ceramic substrates with the conventional three-electrode electrochemical configuration composed of a WE, a CE, and a reference electrode (RE) designed, as shown in Fig. 1.The WE and CE were made of gold, while the RE was made of silver.The WE disk had a diameter of 1.6 mm with all electrodes having a thickness of 100 µm.Two different types of devices with the same geometrical structure were tested: DRP-C223AT (AT) and DRP-C223BT (BT), being the main difference in the temperature of annealing of the conductive gold ink during the fabrication process.The BT devices presented a highly rough surface due to low-temperature ink curing [20].Meanwhile, the AT electrodes had a smoothened surface since they were cured with a high-temperature process [20], [21].

B. Experimental Setup and Measurements
The electrical responses of the sensors were characterized through CV and EIS.Experimental measurements were performed by depositing on the electrodes' surface a 100-µL drop of a phosphate-buffered saline (PBS) solution containing ferro/ferricyanide Fe(CN) 3−  6 /Fe(CN) 4−   6   (FeCN) at different concentrations.The CV measurements were carried out by the CH-440 A potentiostat, sweeping the potential between −0.2 and +0.45 V with scan rates ranging from 30 to 100 mV/s.The EIS measurements were performed in the three-electrode configuration by using the combination of an electrochemical interface, Solartron 1287, and an impedance spectroscopy analyzer, Solartron 1260.The dc bias was set to the equilibrium potential E 0 , which was extrapolated as the middle potential between the cathodic and anodic peaks in the CV curves, and V ac was set to 10 mV with impedance measurements in the frequency span of 1 Hz-100 kHz.EIS measurements with a three-electrode configuration are fundamental to characterize the single WE, avoiding CE contributions, which may influence the results [22], [23].Each test was repeated at least four times on different devices, and the obtained curves were used for the model calibration and validation.

III. PROPOSED MULTIPHYSICS MODEL A. Three-Dimensional Geometry and Meshing
COMSOL Multiphysics software allowed to define finite-element models by combining the equations characterizing different physical processes [24].The COMSOL simulations were based on a 3-D geometry, where each structural component was characterized through specific materials' proprieties defined by the user.
The shapes of each 3-D element of the sensor layout were obtained by using a combination of CAD operations and geometric primitives.The geometrical structure of the devices comprised the ceramic substrate, the three electrodes (WE, CE, and RE), the metal contact tracks, and an insulating layer.The solution drop was geometrically designed as a semi-sphere with a 4 mm radius in contact with the active area of the electrodes.The dimensions and the thickness of each geometrical component were defined as close as possible to the data retrieved from the manufacturer datasheets.The roughness of the surface was not implemented in the 3-D geometrical layout of the model, as geometrical features at microscopic level are not implemented in the geometry.However, the calibration is intended to compensate the inadequacy of the planar electrodes' layout.The 3-D geometries for both the considered devices are reported in Fig. 2(a).
The 3-D structures were finely meshed with tetrahedral elements in all the domains of interest, as shown in Fig. 2(b).Meshing operations were performed on the electrodes, the contact tracks, and the solution drop.Further meshing refinements were required on the solution drop around the electrode boundaries, in order to accurately model the redox reactions and diffusion phenomena occurring at the metal/solution interface.The mesh was further refined in the interface region around the WE, as emphasized in Fig. 2(c).
The first boundary layer was defined with a stretching factor of 1.1 and thickness of 1 µm, which is calculated from where D is the diffusion coefficient, and t is the diffusive minimum time [25], [26].The other layers' thickness proportionally decreases until reaching the electrode surface.
The final mesh had a total of ∼194.000elements, while the computational time of a single simulation was around 1 h for an Intel i9-11900K processor, reaching a satisfying tradeoff between the simulation accuracy and the computational burden.

B. Multiphysics
The simulation comprised a combination of four physics modules.

4) Electrical circuits (cir).
The secondary current distribution (cd) was the core physic, which modeled the electrochemical potential and current distribution in the electrodes, at the metal/solution interface, and in the solution.The conductivity of the supporting electrolyte, i.e., PBS, was set to 1.6 S/m, obtained from the PBS supplier.The electrolyte at the metal/solution interface was modeled considering the double layer capacitance (C dl ), while the electrode/solution interface kinetics was based on the following extended Butler-Volmer equation: where i 0 is the exchange current density, C R and C O are the redox species concentrations, which are the dependent variables of the model calculated at each time point in the mesh network, and α a and α c are the electron direct transfer coefficients, which indicate the symmetry of the redox reaction and represent the quantity of electrons directly transferred across the interface.The overpotential η is defined as s − l − E eq , where s is the electrode potential, l is the electrolyte potential, and E eq is the equilibrium potential of the electrode metal/solution interface.The transport of diluted species (tds) modeled the transport and the reactions of the reduced species and oxidized species in the supporting electrolyte.For each redox species, the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
module implemented the Nernst-Plank equation where D and C are, respectively, the diffusion coefficient and the concentration of the species, z is the charge state, φ is the electric potential in the solution, and v is the flow velocity of the bulk motion of the solution.The last term accounts for convection, which is negligible in our simulations as the solution is static.
The electrons (e − ) flow between electrode/solution interfaces is guided by the redox reaction Ox + e − Red, and the reaction induces the changes in the concentration of the redox species C R and C O at the metal/solution interface depending on the reaction direction.Hence, the tds module was combined with cd by coupling ( 3) and ( 1) to consider the effect of concentration on the electron flux across the interface.
The electric current module (ec) simulated the voltage and current distribution on the conductive tracks and terminals by solving the current conservation equation based on the classical electromagnetism of Maxwell's laws.
The device terminals of ec physic were coupled with the electric circuit module (cir).By the cir module, the electrical circuit of Fig. 3(a) was designed to simulate the behavior of an ideal potentiostat, which is connected to the sensor terminals to generate, condition, and read the electric signals.
The instrument keeps a stable and constant voltage between the WE and RE by regulating the current flowing through the CE, thus avoiding the reference polarization.While works in literature often neglected to consider this element [14], the inclusion of the potentiostat gives a key contribution to the final modeling of the sensors' electrical response, as well as flexibility on simulating different measurement methods.The circuit comprised two ideal operational amplifiers (op-amp), defined as reported in Fig. 3(a), with an input resistance (R in ) of 10 M , an output resistance (R out ) of 50 M , and a voltage-controlled voltage source with a gain of 10 5 (A op ).In the potentiostat circuit, a linear sweep tension (Vlsg) was generated and the current flowed through R1 (1 k as R2) as Vlsg/R1.Since the op-amp A absorbed negligible current at the input terminals and given the op-amp B configuration, the output voltage of B was −Vlsg and the reference current was approximately null.Lastly, the R_OUT was connected to mass and equal to 1 , so the WE was approximately kept at 0 V. Hence, Vlsg controlled the potential between WE and RE, which was kept fixed by the changes in CE voltage.The WE current, flowing through R_OUT, was the output signal of interest, i.e., the current flowing in/out the WE terminal.The applied voltage Vlsg and the current on R_OUT simulated by the potentiostat model are reported in Fig. 3(b).

IV. MODEL CALIBRATION
The proposed model was calibrated through the definition of key parameters so to match the experimental results.The key parameters of the model were the redox species diffusion coefficients D O and D R , the charge transfer coefficients α a and α c , the redox equilibrium potential E 0 , and the exchange charge density current i 0 , which are related to (2), as well as an estimation of C dl .We extrapolated experimentally the diffusion coefficients of ferro/ferricyanide D O and D R from the EIS response at low frequency.In fact, for reversible or quasi-reversible reactions, the Warburg coefficient A W is defined as The value of A W can be retrieved from experimental EIS data fitting with Randles' circuit.Hence, D O and D R were calculated from A W , knowing the bulk concentration of the redox species (C b O and C b R ) and the WE area (A).We compared the obtained diffusivity coefficients with ferro/ferricyanide D O and D R from [22].The experimental values were consistent with their theoretical counterparts; thus, both D O and D R were imposed with the value 7.6 × 10 −6 cm 2 /s.
Both the transfer coefficients α a and α c were defined as 0.5 since in the experimental measurements, the cathodic current peak was comparable with the anodic peak absolute value.E 0 and C dl were retrieved, respectively, from CV and EIS measurements in a solution of 10-mM Fe(CN) 3−  6 /Fe(CN) 4−   6   in PBS (Fig. 4).The CV measurements were carried out with a sweep rate of 100 mV/s and E 0 was calculated Lastly, i 0 was calculated considering only its dependency on the redox species concentration and assuming to work at redox reaction equilibrium potential.The parameter was calculated as i 0 = 5.7 µA/mm 2 for a bulk concentration of 10 mM, i 0 = 2.8 µA/mm 2 for 5 mM, and i 0 = 1.4 µA/mm 2 for 2.5 mM.The outcome showed a good agreement between the experimental data and simulated signals.
Since the simulated data allow the investigation of the electrochemical phenomena on the electrodes, we studied the current density distribution in the electrolyte.As shown in Fig. 5, most of the current density in the electrolyte is distributed between the WE and the CE electrodes and the concentration is strongly dependent on the position and time, which is the expected behavior.More detailed analysis of the uniformity of the spatial concentrations will be the future object of study, aiming to address also the parameters influence on the modeling of electrochemical phenomena.Overall, one of the greatest advantages of the simulation is the investigation of the microscopical electrochemical phenomena at the electrode/solution interface, as spatial and/or temporal heterogeneity in the current densities, concentrations, and reactions of the electrochemical cell, which can be useful to develop optimized sensor configurations.

V. MODEL VALIDATION A. Influence of the Scan Rate
Using the calibrated model, we simulated CV responses at different scan rates.The results were then compared with experimental measurements on AT and BT devices performed at the same scan rates.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.The simulations were performed using the scan rates of 30, 70, and 100 mV/s at a concentration of 10 mM of FeCN in PBS.All the other model parameters were kept fixed.Fig. 6 shows the comparison between the simulated results and the experimental data, reported as mean and standard deviation (m ± std), at 100 mV/s of scan rate for AT and BT devices.Notably, the simulations curves fit well the experimental data.In AT sensors [Fig.6(a)], the simulated cathodic peak current is slightly higher than the experimental peak.This is probably due to a little misalignment of the symmetry coefficients α a and −α c .Meanwhile, the anodic and cathodic peaks of BT devices [Fig.6(b)] greatly overlap the experimental data.
The simulated peak current (I p ) with respect to the experimental m ± std value at different scan rates for AT and BT devices is reported, respectively, in Fig. 7(a) and (b).It is important to highlight that all the simulated I p values follow a linear behavior for both AT (R 2 = 0.99) and BT sensors (R 2 = 0.93).This linearity is proportional to the square root of the sweep rate, proving the consistency with the Randles-Ševčík equation.The AT-simulated anodic I p at 100 mV/s [Fig.7(a)] has very little difference from its experimental counterpart, showing only a 5% decrease with respect to the m ± std value.On the other hand, the cathodic I p reports a shift of 9.1 µA between the experimental and simulated peaks.An analogous behavior can be seen at 70 and 30 mV/s, where the cathodic peak increase is, respectively, 7.3 and 3.1 µA.Meanwhile, at 100 mV/s, the BT anodic I p variation [Fig.7(b)] is 1%, while the cathodic peak increase is less than 5%.A slightly more evident difference can be seen at 70 mV/s since both the simulated anodic and cathodic I p vary by 8 µA from their experimental values.At 30 mV/s, this variation is of 5 µA.This is probably caused by some superficial impurities in the electrodes, which reduce the overall current of the experimental system.Fig. 7(c) shows the simulated anodic and cathodic peaks' potential compared to the experimental m ± std in AT and BT devices.All the simulated peak potentials are almost within the experimental error of the measured data, with variations of a few millivolts from the average values.For example, the simulated AT anodic potential at 100 mV/s shows a slight shift of 25 mV, while the cathodic potential at 100 mV/s differs by 45 mV.The final E of the simulated data is 140 mV, which is comparable to the experimental E = 117 mV.Overall, the fit of the simulated curves on the experimental data was considered qualitatively satisfying for both the considered devices, leaving room for further improvement achievable in future works.Hence, we concluded that the proposed model was capable of simulating CV results at different scan rates with an adequately good approximation.

B. Influence of the Bulk Concentration
We simulated CV signals with different bulk redox concentrations in PBS using the calibrated model.Fig. 8(c) shows the comparison between the simulated peaks' positions and the experimental m ± std in AT and BT sensors.
Almost all the simulated peak potentials show only slight variations, with values within the experimental error.For instance, the simulated cathodic potentials at 10 mM differ less than 20 mV from the experimental average values, and anodic potentials report a maximum difference of 33.8 mV.Interestingly, the simulated peak potentials follow closely the trend of their experimental counterparts varying the redox concentration.In particular, the cathodic potentials retrieved at 5 mM are far lower than the values at 2.5 and 10 mM,  suggesting that this behavior may be characteristic for the redox reaction and devices under test.
Hence, the proposed model proved to be able to simulate CV measurements at different redox concentrations, providing well-fit results with respect to the experimental curves for AT and BT sensors, nicely following their electrochemical characteristics.

VI. CONCLUSION
With this work, we presented the modeling of the CV responses of an electrochemical sensor through multiphysics simulations.The proposed model included several physics describing the electrical, chemical, and electrochemical phenomena characterizing the CV signal of the sensors.The model key parameters were calibrated through dedicated experimental CV and EIS measurements with a solution of 10-mM FeCN in PBS as an electrolyte.The calibrated model simulations were then compared with new experimental CV measurements performed on two commercial sensing devices with different surface roughness.The simulated CV curves fit well the experimental responses even when changing several experimental conditions, such as different redox concentrations and scan rates.The simulated current peaks had values comparable with the experimental error with worst case variation of 37% from the measured values.The model also proved to be useful to investigate the microscopic processes involved in complex electrochemical systems.The approach may be useful to investigate the system limitations, and heterogeneity, and optimize the devices' layout and design.
Future work will aim to test other electrodes with different layouts and proprieties in order to simulate their electrochemical characteristics and optimize the fabrication process, introducing also the modeling of electrodes' surface modification.

Manuscript received 24
May 2024; accepted 1 June 2024.Date of publication 11 June 2024; date of current version 1 August 2024.This work was supported in part by the Catalan Institute of Nanoscience and Nanotechnology (ICN2) through the Catalan Institution for Research and Advanced Studies (CERCA) Program and the Generalitat de Catalunya, in part by ICN2 through the Severo Ochoa Centres of Excellence Program under Grant SEV-2017-0706, in part by the Severo Ochoa Centres of Excellence Program under Grant CEX2021-001214-S and Grant MCIN/AEI/10.13039.501100011033.The associate editor coordinating the review of this article and approving it for publication was Dr. Satyendra K. Mishra.(Corresponding author: Stefano Bonaldo.)Stefano Bonaldo, Lara Franchin, and Alessandro Paccagnella are with the Department of Information Engineering, University of Padova, 35131 Padua, Italy (e-mail: stefano.bonaldo@unipd.it).

Fig. 1 .
Fig. 1.(a) Layout of the devices under test.(b) Drop of redox mediator FeCN in electrolyte solution on the electrodes under test to perform the experimental measurements.

Fig. 2 .
Fig. 2. (a) General view of the 3-D geometrical structure of the devices under test in COMSOL Multiphysics.(b) Mesh of the model.(c) Magnified view of the mesh refinement at the electrode/solution interface around the WE.

Fig. 3 .
Fig. 3. (a) Schematic of the potentiostat electrical circuit and the op-amp circuit (inset).(b) Triangular voltage signal and the current signal simulated by the potentiostat electrical circuit.

Fig. 4 .
Fig. 4. Experimental measurements of the CV and EIS responses of the devices with 10-mM FeCN in PBS.The CV measurements are carried out to retrieve the equilibrium potential E 0_W and E 0_C of the WE and CE electrodes.(a) and (d) CV responses at scan rate of 100 mV/s using the WE electrode of AT and BT devices.(b) and (e) CV responses at a scan rate of 100 mV/s inverting the CE and WE electrode of AT and BT devices.(c) and (f) Three-electrodes EIS curves of AT and BT devices.Dashed lines indicate the experimental measurements, while continuous lines indicate the data fit.The equivalent electrical circuit used for the data fit of EIS measurements is shown in inset.

Fig. 5 .
Fig. 5. Simulation of the current density distribution on the electrodes' surface (gradient-colored areas) and the current direction (black arrows) in the solution at the anodic peak of the CV response.

Fig. 7 .
Fig. 7. Simulated CV current peaks (red) compared with mean and standard deviation (blue) of experimental current peaks at different scan rates in 10 mM Fe(CN) 3− 6 /Fe(CN) 4− 6 in PBS.The black line shows the linear fitting.(a) AT sensors.(b) BT sensors.(c) Comparison between simulated and experimental CV anodic (top) and cathodic (bottom) peak potentials of AT and BT devices at different scan rates.

Fig. 8 .
Fig. 8. Simulated CV current peaks (red) compared with mean and standard deviation (blue) of experimental current peaks at 100-mV/s scan rate at different Fe(CN) 3− 6 /Fe(CN) 4− 6 concentrations.The black line shows the linear fitting.(a) AT sensors.(b) BT sensors.(c) Comparison between simulated and experimental CV anodic (top) and cathodic (bottom) peak potentials of AT and BT devices at different concentrations.