Numerical modeling of electrochemical transport and discharge mechanisms in hierarchically porous lithium-ion electrodes

Abstract

Hierarchical porous microstructures are promising candidates for lithium battery electrodes that maintain capacity and specific energy at high charging and discharging rates. The electrolyte channels embedded in these structures facilitate Li transport throughout the active material such that the accessible capacity is enhanced. We performed numerical simulations of the discharge process for four modeled electrode structures in order to investigate the impact of the size of electrolyte channels on the electrochemical performance. The results show that the size ratio of the electrolyte channels and active materials columns determines the discharge characteristics of hierarchical porous electrodes. Depending on the size ratio, electrodes discharge from the separator side and the current collector simultaneously which can enhance the accessible capacity and specific energy. The results from our simulations can aid in designing tailored hierarchical porous electrode structures for fabrication of electrodes with enhanced capacity and rate capability.

---

1. Introduction

Li ion batteries are the primary energy storage devices for a wide range of e-technologies such as electric vehicles and mobile devices. With ca. 30% of total U.S. energy consumption used by the transportation sector in 2023,¹ improved battery technology will be critical for transporting people and goods in the future to electrify the transportation sector. However, high costs, limited lifetime, and limited energy density of Li ion electrodes hold back large-scale applications of Li batteries.^2–4 High specific energy is nearly always desirable in battery systems, but it is especially important in batteries for electric vehicles. Current electric vehicle batteries achieve specific energy between 148 and 158 Wh kg⁻¹ with cell-level specific energy ranging from 200 to 235 Wh kg^−15,6 but further advances are needed for larger market penetration, particularly for battery electric vehicles (BEV). For example, the DOE ARPA-E program PROPEL-1K set a specific energy target of 1000 Wh kg⁻¹.⁷

One way to improve specific energy for a particular combination of anode and cathode active material is to maximize the mass fraction of active material, which means minimizing the fraction of inactive material such as the separator, the binder, electrode additives for improving electronic conductivity, and the current collector. A seemingly simple strategy for maximizing the mass fraction of active materials is to make the electrodes as thick as possible.^8,9 However, the thickness of electrodes supported on current collector foils cannot be increased without bound due to limits on transport rates of electrons and lithium ions in the electrodes.^10,11 Because of these limitations, current battery cathodes tend to be limited to areal mass loadings of 20–40 mg cm⁻² with average porosities between 20–40%, which results in electrode thicknesses that are usually less than 80 micrometers.^12–15 Attempts to increase mass loading and electrode thickness beyond these limits without optimizing the electrode microstructure usually result in substantial portions of the active material being inaccessible except at quite low charge/discharge rates that would be unacceptable for vehicle use.

Liquid electrolyte mass transport is a major limiting factor for the performance of high-power Li ion batteries.⁶ At high discharge rates, the slow transport of Li ions leads to Li depletion near the current collector resulting in a drop of cell voltage. Ion transport in thick electrodes can be particularly problematic. A typical slurry-casted battery electrode with a porosity near 30% contains much less electrolyte within the battery pores than is required to fully charge or discharge the battery. Furthermore, the pore structures and morphology are random, leading to higher tortuosity for ion transport.^16,17 Salt diffusion into and through the electrode pore space is slow and can easily limit the battery charge/discharge rate. Optimization of porous electrodes by introducing hierarchical anisotropic porosity is a promising way to overcome the transport limitations and improve battery performance.^18,19

Experimental work has highlighted the utility of directional freeze casting for creating battery electrode microstructures with enhanced transport characteristics.^18,20 In this approach, an electrode slurry is rapidly frozen, creating ice crystals that upon subsequent removal by sublimation leave behind unidirectional columnar macropores that may be filled with electrolyte. Azami-Ghadkolai et al.¹⁸ have combined directional freeze casting with tape casting to develop a process that can be easily adopted to make Li-ion battery electrodes at a commercial scale. Using Mo-doped lithium titanate (MoLTO, Li₄Mo_xTi_5−xO₁₂) electrodes with a range of microstructures using both conventional tape casting and freeze tape casting were fabricated. Cells were fabricated using these electrodes as cathodes vs lithium metal anodes. Galvanostatic charge–discharge testing showed that the freeze-tape-cast samples exhibit higher specific capacity under all conditions but particularly under high-rate charging/discharging. These results likely reflect the greater access of electrolyte to the battery active material in these electrodes compared with conventional tape-cast electrodes with the same mass loading.¹⁸ A deeper understanding of the dependence of ion transport characteristics on the materials properties (electric conductivity, ionic conductivity) and microstructure (porosity, tortuosity, anisotropy) of Li ion electrodes will be instrumental in accelerating the development process for engineered porosity electrodes.

Theoretical models of Li ion batteries originate from the pioneering work of Newman and co-workers using concentrated electrolyte theory.^21–24 Various continuum theories have been put forward in the literature.^8,25–31 A popular approach is the 1 + 1D model,⁸ a pseudo-2D model at the cell-level that is derived by using systematic volume averaging over the porous microstructure. The transport in the electrolyte is governed by a macroscopic mass balance that describes diffusion and migration of Li ions. The migration flux is driven by a potential gradient that is obtained from a charge balance. The transport in the active material is described by diffusion of lithium. The de-/intercalation of Li ions is represented by a source term in the mass balance weighted by the specific surface area. The 1 + 1D model is computationally efficient and can be solved on a standard desktop computer. While the pseudo-2D model can be used for simple geometries, it is not directly applicable to hierarchically structured porous electrode materials.³²

A multiscale approach is required to model the charge transport in hierarchical porous electrodes. While multiphase porous electrode theory has been developed to capture phase separation and intraparticle thermodynamics in active materials,³³ the present work focuses on transport limitations arising from hierarchical electrode geometry and anisotropic electrolyte pathways. The overall charge/discharge can be modeled at the macroscale level, while the spatial distribution of Li ions in the microstructure is either modeled explicitly,³⁴ or by means of a surrogate model such as a subgrid model for Li diffusion within active particles. Unlike transmission line models that describe porous electrodes using lumped or distributed equivalent circuits,³⁵ the present work employs continuum transport equations that explicitly resolve concentration and potential fields in hierarchical electrode geometries. In this context, Kashkooli et al.³⁶ have investigated the effect of particle size of monodisperse active particles on the performance of LTO electrodes. They demonstrated that a Newman pseudo-2D model can be employed to simulate various electrode materials when the model parameters are set to experimentally measured values of the actual electrode material.³⁶ Rashid et al.³⁷ presented numerical simulations of LTO half-cells and investigated the effect of particle size, lithium diffusivity, and electrode thickness on the capacity of Li-LTO cells. Their results indicated that larger particle size and lower lithium diffusivity limit the utilization of active material, leading to lower specific capacity at higher discharge currents. Allu et al.³¹ employed a volume averaged formulation of the transport equation to model 3D prototype electrodes. They showed that the volume averaged approach is capable of capturing spatio-temporal variations of potential and lithium distribution in interdigitated electrode geometries. Cobb and Bianco³⁸ conducted simulations of 2D cross-sections of LiCoO₂ electrodes fabricated by co-extrusion printing. They investigated the impact of the thickness and aspect ratio of the co-extruded domains on the electrode performance and demonstrated that a careful design of the geometry enhanced material utilization in thick electrodes on the order of 150–300 µm. Azami-Ghadkolai et al.³² employed a similar design principle to MoLTO electrodes. By comparing experimental measurements for freeze-tape cast anisotropic porous electrodes with numerical computations for a pseudo-2D model for hierarchical electrode structures, they demonstrated that enhanced material utilization and capacity is facilitated by Li salt transport through electrolyte channels. Their simulation results suggest that freeze tape cast electrodes of thickness up to 750 µm may be fully discharged at a rate of 1C, whereas at the same mass loading, normal tape cast electrodes without electrolyte channels can only be discharged up to 300 µm thickness.

In this work, we investigate the impact of the channel width and aspect ratio on the specific capacity and power density of thick freeze-tape case MoLTO electrodes. We employed numerical simulations of the volume averaged equations corresponding to Newman's model to 2D cross-sections of columnar electrode geometries, using the experimental measurements of Azami-Ghadkolai et al.^18,32 to set the materials parameters. We present results of electrode discharge simulations at varying C rates for three different columnar geometries with the same mass loading. We analyze the dependence of specific capacity and power density on the hierarchical distribution of active material. Our results provide insights into the factors that influence the performance of thick MoLTO electrodes. We discuss how these insights can be leveraged to infer design rules for optimization of the hierarchical porous electrode microstructure.

2. Methods

2.1. Electrokinetic transport in porous electrodes

Theoretical models for charge transport in porous electrodes date back to the work of Newman and coworkers.^21,22,39,40 Newman's model considers the charge transport in the liquid and solid phases by applying concentrated solution theory to porous electrodes. Recently, Lai and Ciucci⁴¹ have revisited Newman's model and re-derived the transport equations via formal volume averaging of the generalized Poisson–Nernst–Planck (PNP) equations. This approach provides a systematic upscaling procedure of the microscopic transport equations to a coarse-grained representation of the porous microstructure at the electrode scale. The volume-averaged equations make it possible to obtain numerical solutions for the charge and discharge characteristics of hierarchically porous electrodes. In this work, we have used the Battery Design Module of the COMSOL Multiphysics software to perform finite element simulations of galvanostatic electrode discharge. In this section, we briefly review the derivation of the coarse-grained transport equations following the presentation in ref. 41.

2.1.1. Poisson–Nernst–Planck equations. At the pore scale of an electrochemical system (but still in the continuum limit), the transport of ionic species is described by the Poisson–Nernst–Planck equationswhere c_k and i⃑_k are the concentration and current density of species k, and ϕ is the electric potential. The current densities i⃑_k = −σ_k∇φ_k are driven by the electrochemical potential where μ_k is the chemical potential, whose concentration dependence is expressed in terms of the activity f_k(c_k)The current densities can then be expressed in terms of the chemical diffusivity as

i⃑_k = −z_keD_k∇c_k − σ_k∇ϕ. (4)

For a monovalent electroneutral salt, the concentrations of cations and anions are equal (c₊ = c₋ = c), such that the current densities become

i⃑₊ = −σ₊∇φ₊ = −eD₊∇c − σ₊∇ϕ, (5)

i⃑₋ = −σ₋∇φ₋ = eD₋∇c − σ₋∇ϕ. (6)

The electric potential can be eliminated from the current densities to express them in terms of the cation electrochemical potential

i⃑₊ = −σ₊∇φ₊ = −t₊σ∇φ₊, (7)

where we have introduced the common definitions of the total conductivity σ, ambipolar chemical diffusivity D, and cation transference number t₊:

σ = σ₊ + σ₋, (9)

D = t₊D₋ + (1 − t₊)D₊, (10)

The total current density thus becomesThis expression retains the cation electrochemical potential gradient and concentration gradient as the driving forces of the charge transport. In the electroneutral formulation, the Poisson equation is replaced by the condition ∇·(i⃑₊ + i⃑₋) = 0.

2.1.2. Charge transport in electrolyte solution. Only two of the current densities i⃑₊, i⃑₋, and i⃑ are independent and it is convenient to replace the anion current i⃑₋ by the total current i⃑. The cation flux and current in the electrolyte solution are written asand the mass and charge conservation equations for the liquid electrolyte phase are

∇·i⃑ = 0. (16)

In the above equations, the subscript l refers to the liquid phase.

2.1.3. Charge transport in active material. The solid electrode phase of the electrode is composed of active particles, where the charge carriers in the solid phase are Lithium ions and electrons. To maintain charge neutrality, the total current density in the active particles must be zero. The flux of lithium ions in the solid phase is thus given by

j⃑_s = −D∇c_s. (17)

Furthermore, it is assumed that in the current collector the electron concentration is sufficiently high such that ∇·i⃑_s = 0. The mass and charge conservation equations then become

∇·(σ_s∇φ_s) = 0. (19)

2.1.4. Lithium intercalation. The intercalation of lithium is driven by the difference in the electrochemical potentials between the liquid and the current collector and the equilibrium potential φ₊ + φ₋ − φ_eq. The intercalation current density at the surface of active particles can be written in the formThe prefactors i^l–s and i^s–l are proportional to the concentration of lithium in the liquid and solid that can contribute to the exchange reaction

i^s–l = c_si^s–l₀, (22)

where c_s,0 indicates the maximum concentration of intercalation sites in the solid and i^l–s₀ and i^s–l₀ are constants independent of concentration. The exchange current density can be recast in the formwhereeqn (23) is commonly written in the form of the Butler–Volmer equation³⁵where α_c and α_a are symmetry factors, V_OC is the open circuit voltage, and i₀ is a constant. It is worth noting that the comparison of (23) and (26) involves a slightly different interpretation of the potentials.⁴¹ For a conceptual discussion of the meaning of the open-circuit voltage, we refer the reader to ref. 42.

2.2. Volume averaged transport equations for porous electrodes

To model hierarchical porous electrodes, the microscopic porosity is treated as an effective porous medium. It is assumed that the microscopic features of the electrode material are much smaller than the length scale of interest such that its effective properties can be obtained by averaging over a representative volume element (REV).^39,41 The volume average of a quantity f over a volume with porosity is defined as²¹Averages of the time derivative, gradient, or divergence can be obtained in a similar fashion and turned into derivatives of the volume averaged quantities with surface integral corrections. Integration along a tortuous path dz = τdx with the tortuosity τ yields the volume-averaged quantities⁴³

2.2.1. Liquid electrolyte phase. In the liquid electrolyte, the concentration, flux, and current are defined within the liquid volume fraction . Volume averaging the flux (13) and the current (14), we obtainThe surface integral turns into a source term that represents the exchange of lithium between the liquid electrolyte and the solid electrode. It only enters the cation flux since there is no exchange of anions. The volume-averaged divergence terms arewhere and a_v = A/V is the specific surface area of electrode active material. The latter can be expressed in terms of the radius r_p and volume fraction of active particles as . The volume averaged mass and charge conservation equations in the electrolyte phase can then be writtenwith the average fluxes given in (31) and (32).

2.2.2. Solid electrode phase. Different approaches for treating the solid electrolyte phase have been proposed in the literature. In principle, volume averaging can be applied to the solid phase as was discussed by Allu et al.³¹ This leads to the volume-averaged equationswith the lithium flux and electronic currentHowever, the diffusion of lithium in the active material is much slower than in the liquid electrolyte such that the volume averaging is inaccurate due to the inhomogeneous concentration. One way to retain the concentration gradients in the active material in the model is to keep the microscopic equations and explicitly solve them on a finer scale, e.g., on a subgrid.

∇i⃑_s = −i_loc, (42)

with the lithium flux and electric currentThis gives rise to a multiscale model that couples the microscopic evolution of the solid phase concentration c_s to the mesoscopic evolution of the volume-averaged electrolyte concentration 〈c_l〉^l. The coupling is represented by the exchange current i_ex, which is explicitly calculated with the surface concentration in the active particles and enters the electrolyte equations as the surface integral i_loc.

2.3. System of porous electrode equations

The system of coupled micro-macro volume averaged equations is obtained by collecting the equations in the following form,⁴¹ where we have dropped the brackets for the intrinsic volume average in the electrolyte phase.

∇·i⃑_l = i_loc, (47)

∇·i⃑_s = i_loc. (48)

The fluxes and currents are

j⃑_s = −D^eff_s∇c_s, (50)

i⃑_s = −σ^eff_s∇φ_s, (52)

and the exchange current density is given byThe effective transport coefficients in eqn (49)–(52) are

2.4. Boundary conditions

A typical electrode system consists of the electrodes, a separator, and the current collectors. In this work, we restrict ourselves to a cathode that is connected to a separator on one side and bounded by a current collector on the other.^18,21,39 The second boundary of the separator is bounded by an anode surface. To model the components, boundary conditions have to be applied at the interface between the separator and the cathode, and at the current collector and the anode surface. The boundary conditions are expressed in terms of the fluxes/currents and, for the purposes of computation, rewritten in terms of the derivatives of concentration and potential.

2.4.1. Interfaces between porous electrode and separator/electrolyte. At the surfaces of the active material, lithium ions can cross the electrode domain boundary through the electrolyte and the conservation laws require continuity of j⃑_l and i⃑_s. The current in the solid phase is subject to an insulation boundary condition i⃑_s = 0.

j⃑_{l,active material} = j⃑_l,separator, (55)

i⃑_{l,active material} = i⃑_l,separator, (56)

i⃑_s = 0. (57)

2.4.2. Current collector. The current collector is a solid conductor and hence the electrolyte fluxes have to vanish. The applied external current i_app imposes a boundary condition on the current in the solid phase.

n̂·j⃑_l = 0, (58)

n̂·i⃑_l = 0, (59)

n̂·i⃑_s = −i_app. (60)

2.4.3. Electrode surface reaction. Instead of an explicit anode domain, an electrode surface reaction is used to model the intercalation of lithium ions in the anode materials. The exchange reaction imposes a surface flux i_electrode that is given in the form^39,41where and . The lithium flux and total ionic current in the electrolyte then are

n̂·i⃑_l = i_elec. (63)

3. Modeling hierarchical anisotropic porous electrodes in COMSOL

In this work, we use the Battery Design Module of the COMSOL Multiphysics software to perform finite element simulations of galvanostatic electrode discharge. The columnar structure of freeze-tape cast electrodes is modeled as an array of rectangular domains of microporous active materials interdigitated by rectangular electrolyte channels. The 2D cross section of such a hierarchical electrode structure is illustrated in Fig. 1. We varied the size of the electrolyte channels and the active materials columns while keeping the overall width of the electrode fixed at W_tot = 480 µm. Specifically, we created geometries with the size ratios of the electrolyte channel and active material w_elec : w_LTO set to 10 µm : 10 µm (equisized), 10 µm : 20 µm (narrow channels), and 20 µm : 10 µm (wide channels). To keep the mass loading fixed at 30.329 mg cm⁻², we assumed a constant electrode area of 70.84 mm² and adjusted the thickness h_LTO of the electrode (455 µm for equisized channels, 340.88 µm for narrow channels, 681.75 µm for wide channels). The geometric parameters of the model are summarized in Table 1.

Fig. 1 Schematic of the hierarchical anisotropic porous electrode structures. The model assumes regular aligned active materials columns of width W_LTO that are separated by electrolyte channels of width W_elec. The total width of the electrode was kept at W_total = 480 µm. The thickness of the active material is denoted by h_LTO.

Table 1 The geometric parameters of the 2D cross-section used to model varying width of the electrolyte channels while keeping the mass loading fixed

Parameter	Variable	Value
		No channels	Narrow channels	Equisized channels	Wide channels
# of MoLTO columns	n	1	16	24	16
Column width	W LTO	480 µm	10 µm	10 µm	20 µm
Channel width	W elec	0 µm	20 µm	10 µm	10 µm
Electrode width	W tot = n(WLTO+ Welec)	480 µm
Electrode thickness	h LTO	250 µm	340.88 µm	454.5 µm	681.75 µm
Separator thickness	h sep	25 µm
Out-of-plane area	A	70.84 mm^2

---

We consider freeze-tape cast electrodes made of MoLTO (Li₄Mo_0.3Ti_4.7O₁₂) as the active material and 1 M lithium hexafluorophosphate (LiPF₆ 3 : 7 EC : EMC) as the electrolyte. A battery half-cell consists of the electrode sandwiched between a polymer separator (25 µm microporous PP/PE/PP membrane) and the current collector. The other surface of the separator is attached to Li metal as the opposite electrode. The materials properties of these components are based on experimental measurements³² and are listed in Table 2.

Table 2 Materials parameters and electrochemical cell parameters were set based on experimental measurements by Azami-Ghadkolai et al.²⁷ to match the properties of freeze tape-case MoLTO electrodes

Property	Variable	Value
Active material (MoLTO)
Diffusion coefficient	D s	1× 10^−17 m^2 s^−1
Conductivity	σ s	1 S m^−1
Equilibrium potential	E eq	Interpolated
Volume fraction of active material in each column	ε s	0.3891
Porosity of active material	ε l,LTO	0.4768
Tortuosity of active material	τ l,LTO	5
Density of active material	ρ LTO	3.43 g cm^−3
Grain particle radius	R p	100 nm
Theoretical capacity	Q s	169.8 mAh g^−1
Maximum SOC	SOCmax	0.944
Minimum SOC	SOCmin	0.003
Rate constant of active material	θ	5 × 10^−11 m s^−1
Reference exchange current	i 0,ref	0.098964 A m^−2
Electrolyte (LiPF6 3 : 7 EC : EMC)
Diffusion coefficient	D l	Interpolated
Conductivity	σ l	Interpolated
Transference number	t +	Interpolated
Activity coefficient	f cl	Interpolated
Separator
Porosity of separator	ε l,sep	0.39
Tortuosity of separator	τ l,sep	4.02
Lithium metal
Reference exchange current	i 0,surf	1.2482 A m^−2

---

To simulate a galvanostatic discharge cycle, the half-cell is subjected to an applied current at the current collector corresponding to the prescribed C rate, while the potential of the Li metal surface is kept constant at ϕ_s,ext = 1 V. The initial Li concentrations are set to c_l,init = 1000 mol m⁻³ and c_s,init = 68.96 mol m⁻³ (corresponding to the maximum state of charge of MoLTO). The electrochemical cell parameters are listed in Table 3.

Table 3 Electrochemical parameters for simulations of electrode discharge. The values are based on experiments reported by Azami-Ghadkolai et al.²⁷

Property	Variable	Value
Active material (MoLTO)
Initial lithium concentration	c s,init	SOCmax × ρ × Qs/F
Initial potential	ϕ s,init	0 V
Electrolyte (LiPF6 3 : 7 EC : EMC)
Initial concentration	c l,init	1000 mol m^−3
Initial potential	ϕ l,init	−2.6 V
Lithium metal
External electric potential	ϕ s,ext	0 V

---

The numerical integration of the governing equations was performed in COMSOL using the PARDISO solver. All geometries were meshed using square grid cells of size 2 µm. The maximum time step of the time-dependent solver was limited to 10 s and the solver tolerance was set to 10⁻⁴. This value was chosen based on preliminary simulations probing stability and accuracy at 2 µm spatial resolution. The discretization parameters used for the numerical solver in COMSOL are listed in Table 4. The simulations were conducted with COMSOL 5.6 on Clemson University's Palmetto cluster using 56 cores of an Intel Xeon 6238R CPU and 300 GB of RAM.

Table 4 Discretization parameters used to tune the numerical solver in COMSOL

Parameter	Variable	Value
Mesh size	dx	2 µm
Output interval	dt	500
Maximum time step	maxstep	10 s
Tolerance	tol	1 × 10^−4

---

4. Results and discussion

We simulated the discharge process of four different electrode structures: a normal tape-cast electrode without electrolyte channels, and three freeze tape-cast electrodes with three size ratio of liquid electrolyte channel: active electrode wall thickness: narrow (w_elec : w_LTO = 10 µm : 20 µm); equisized (w_elec : w_LTO = 10 µm : 10 µm), and wide (w_elec : w_LTO = 20 µm : 10 µm). All simulations were run until the cell potential decayed to a cutoff voltage of 1 V, and the time to reach this condition was recorded as the discharge time. The rate of discharge is commonly specified as the C rate, where 1C corresponds to the discharge current at which the maximum capacity would be discharged in one hour. Based on the experimentally measured properties given in Table 2, we calculated the applied current density for a prescribed C rate relative to the theoretical maximum capacity as

4.1. Electrode discharge

Fig. 2 shows the evolution of the cell potential over time for electrodes with a mass loading of 30.329 mg cm⁻² and at discharge rates of 1C and 2C. At the lower discharge rate of 1C, all four electrode structures can be fully discharged in the theoretical discharge time of 1 hour. The equisized and narrow channel geometries reach the cutoff voltage of 1 V at similar times, while the discharge stops earlier for the geometries with wide channels and with no channels. This indicates that the equisized and narrow channels electrodes provide nearly the theoretical capacity, while the electrodes with wide channels and without channels suffer from some loss of capacity. This effect becomes more pronounced at the higher discharge rate of 2C, where the electrodes with equisized channels and narrow channels can still be discharged for close to 30 minutes, while the electrodes with wide channels and without channels shut off earlier. Interestingly, at 2C, the electrode structure with narrow channels provides a lower capacity than the electrode without channels. We will investigate this observation in more detail in Section 4.2.

Fig. 2 Development of the electrode half-cell potential over time during discharge at rates of 1C and 2C for the four electrode structures without channels and with narrow (10 µm : 20 µm channel to column width ratio), equisized (10 µm : 10 µm channel to column width ratio) and wide channels (20 µm : 10 µm channel to column width ratio).

Our results are in line with the experimental observations of Azami-Ghadkolai et al.³² that the electrolyte channels in freeze tape-cast electrodes can enhance the accessible capacity and thus improve the rate capability. It is worth noting, however, that the simulations in ref. 32 have not been performed for the experimentally measured materials properties of MoLTO. This is evident in the concentration plots in Fig. 4 and 5 of ref. 32, where the maximum state of charge (SOC) exceeds the value of SOC_max = 0.9 reported in Table 2 in ref. 32. The accuracy of the previous numerical results is thus an open question. While we here report results based on the actual materials properties of MoLTO, due to the unavailability of experimental data, the comparison to experimental measurements is left for future work.

4.2. Rate capability

To evaluate the rate capability of the four different electrode structures, we varied the discharge rate from 1C to 8C. For each simulation, we recorded the discharge time at which the cell potential dropped to the cutoff voltage, and calculated the areal capacity

The dependence of the areal capacity on C rate is plotted in Fig. 3 for the four electrode structures at the same mass loading of 30.329 mg cm⁻² corresponding to a theoretical capacity of 5.15 mAh cm⁻². The results show that the rate capability of hierarchically anisotropic electrodes can be substantially enhanced due to the presence of electrolyte channels. At low discharge rates up to 1C, all four electrodes show a limiting capacity of just under 5 mAh cm⁻². For the electrodes with wide channels and without channels, the capacity decreases with increasing C rate showing a notable capacity loss at higher discharge rates. The electrode with equisized channels maintains a high capacity up to 2C, while the electrode with narrow channels can be discharged at up to 3C without substantial capacity loss. In comparison, electrodes with equisized channels and narrow channels show a better rate capacity than electrodes with wide channels and without channels. The best rate capacity is observed for electrodes with narrow channels. Fig. 3 also confirms that electrodes with wide channels show an even lower rate capability than the electrodes without channels. This indicates that the size ratio of active material columns and electrolyte channels in freeze tape-case electrodes is critical to the discharge performance and provides a means to optimize the fabrication process.

Fig. 3 Rate capability of the four electrode structures at the same mass loading of 30.329 mg cm⁻². The thickness h_LTO of the electrode was 250 µm without channels, 455 µm for equisized channels, 340.88 µm for narrow channels, and 681.75 µm for wide channels. The superior rate capability of the equisized and narrow channel electrodes can be observed.

The impact of the electrolyte channel width on the accessible capacity can be rationalized by analyzing the concentration profiles of Li ions throughout the electrode geometry. Fig. 4 shows the concentration profiles of salt concentration in the electrolyte and the state of charge of the active material throughout the electrode for the four electrode geometries at 1C. The state of charge is also shown as a 2D plot to illustrate that the major concentration variations occur in the thickness direction. The concentration profiles for the electrolyte show that a salt concentration gradient develops, where the salt concentration is highest at the separator and lowest at the current collector. This gradient is a result of the consumption of Li through intercalation within the microporous active material. To replenish the salt concentration, salt ions have to be transferred from the separator side of the electrode which is limited by the transport properties of the electrolyte. For the electrode without channels, the Li salt near the current collector is increasingly depleted to a near-zero concentration at the end of discharge. The intercalation of Li thus becomes entirely limited by the salt transport in the electrolyte. For the electrodes with channels, the surplus of salt in the electrolyte serves as a reservoir such that the depletion is less pronounced. Azami-Ghadkolai et al.³² have hypothesized that a higher salt concentration decreases the local concentration overpotential such that the Li intercalation rate is greater and the electrode discharges faster in regions of high salt concentration. Therefore, electrodes without electrolyte channels discharge preferentially from the separator side towards the current collector, as indicated in the corresponding SOC profile in Fig. 4. However, the picture becomes different for the electrodes with electrolyte channels. Indeed, the SOC profiles for wide and equisized channels show that the electrodes discharge preferentially from the current collector towards the separator side. This can be attributed to the effect of ohmic losses in the active materials columns. Since the electrolyte channels supply sufficient Li salt, the discharge rate leads to higher current densities and higher ohmic losses. This gives rise to a potential gradient in the active material which lowers the local overpotential near the current collector. Therefore, the electrodes with electrolyte channels can discharge more rapidly at the current collector. At the end of discharge, both regions near the current collector and near the separator have a high SOC, while a region of lower SOC remains inside the active material columns which reduces the accessible capacity.

Fig. 4 Spatial 1D profiles of salt concentration in the electrolyte (left column) and the state of charge of the active material (middle column) throughout the electrode for the four electrode geometries at 1C. The state of charge at the end of discharge is also shown as a 2D plot (right column) to illustrate that the major concentration variations occur in the thickness direction. The state of charge reveals that the location of the inaccessible active material shifts towards the center with decreasing width of the electrolyte channels.

Interestingly, the location of this region changes with the size ratio of the electrolyte channels. For the electrodes with wide channels, the inaccessible region is located near the separator, consistent with the discharge progressing from the current collector towards the separator side. For the electrode with equisized channels, the inaccessible region is located further away from the separator, and for the electrode with narrow channels the inaccessible region is found approximately in the middle of the electrode columns. This indicates that the electrodes with equisized and narrow channels discharge from both sides inwards. This can be attributed to the interplay of the salt concentration gradient in the electrolyte and potential gradients in the active material columns. Furthermore, narrow channels lead to a more frequent alternation of electrolyte and active material domains across the electrode width, which reduces the lateral distance over which Li⁺ must diffuse from the electrolyte reservoir to the active material. As a result, the electrolyte channels act more effectively as distributed Li⁺ reservoirs. Both effects decrease the local overpotential, leading to a greater driving force for discharge in the respective regions. The net effect is thus a competition of discharge near the separator with discharge near the current collector, and the shift of the inaccessible region depends on the relative magnitude of the effect. It can be noted that this effect is associated with the design of the electrodes: to maintain the same areal mass loading, electrodes with wide channels have to be made thicker while electrodes with narrow channels can be thinner. The thickness of the electrodes has a direct impact on the salt transport and the ohmic losses in the active materials columns. The higher ohmic losses along the length of the active material columns lead to a lower discharge rate near the separator, such that the thicker electrodes discharge preferentially at the current collector and the inaccessible region shifts toward the center of the active material columns. This insight can be leveraged to fabricate engineered porosity electrodes with carefully designed channel width and thickness that optimize the accessible capacity.

4.3. Areal capacity: dependence on electrode thickness/mass loading

We further investigated the effect of electrode thickness on the areal capacity and specific energy of the four electrode structures. We varied the thickness of the electrode in the range from 200 µm to 1000 µm and simulated discharge rates from 0.5C to 3C for each thickness. Fig. 5 shows the areal capacity as a function of thickness and mass loading, respectively, for each of the electrode structures. We have used the electrode with equisized electrolyte channels as a reference for the comparison. Fig. 5 clearly shows the maximum thickness at which an electrode can be discharged without loss of capacity. The results show that an electrode without electrolyte channels may be almost fully discharged up to a thickness of approximately 300 µm at 1C, 200 µm at 2C, and 100 µm at 3C, while an electrode with equisized channels can be discharged up to a thickness of approximately 700 µm at 1C, 500 µm at 2C, and 400 µm at 3C. Electrodes with wide channels show a similar trend, however, the areal capacity is lower than for equisized channels. In contrast, electrodes with narrow channels provide a higher capacity compared with equisized channels, but the maximum thickness is slightly smaller (600 µm at 1C, 400 µm at 2C, and 340 µm at 3C). The difference in areal capacity stems from the fact that the width of the active columns is different such that the mass loading of the electrodes with channels increases disproportionately with increasing thickness. To incorporate this effect, we also plotted in Fig. 5 the dependence of the areal capacity on mass loading. The direct comparison shows that electrodes with narrow channels allow the highest mass loading without capacity loss, followed by electrodes with equisized channels. Conversely, electrodes with narrow channels and without channels allow less mass loading than equisized channels. The results indicate that engineered porosity electrodes facilitate higher mass loading than random porosity electrodes. Fig. 5 also shows that it is important to distinguish the independent variable based on which the comparison is made, in this case thickness and mass loading, respectively.

Fig. 5 Comparison of the areal capacity of the four electrode structures as a function of thickness (left column) and mass loading (right column). The equisized channel electrode was used as reference in comparison with the electrodes without channels, with wide channels, and with narrow channels.

4.4. Specific energy: dependence on electrode thickness/mass loading

Another important quantity for assessing the performance of battery electrodes is the specific energy. The specific energy is influenced by both the accessible capacity and the development of the cell potential over time. For our simulations, we obtained the specific energy from integral of the cell potential (cf.Fig. 2) over time
Fig. 6 shows the specific energy as a function of thickness and mass loading, respectively, for each of the three channel porosity structures in comparison with the no channel electrode. The results show that upon increasing the thickness, electrodes with electrolyte channels can deliver a higher specific energy than electrodes without channels. With increasing C rate, the ability of thick electrodes to deliver their maximum specific energy reduces substantially. The fading of specific energy is most pronounced for electrodes with wide channels and without channels. As before, we also plot the specific energy as a function of mass loading in Fig. 6. This comparison clearly shows that at the same mass loading and C rate, electrodes with narrow channels deliver higher specific energy than electrodes with equisized channels, which in turn deliver higher specific energy than electrodes with wide channels or without channels. For the electrodes with wide channels, the ability to deliver high specific energy is reduced due to the larger thickness required to maintain mass loading, which reduces the accessible capacity as discussed in Section 4.2.

Fig. 6 Comparison of the specific energy of the four electrode structures as a function of thickness (left column) and mass loading (right column). The equisized channel electrode was used as reference in comparison with the electrodes without channels, with wide channels, and with narrow channels.

4.5. Normalization with penetration depth

The plots in the previous sections show that the comparison of areal capacity and specific energy depends on whether the results are plotted as a function of thickness or as a function of mass loading. In addition, there can be variations in how the capacity is quantified, e.g., areal capacity vs. volumetric capacity vs. specific capacity. It is thus desirable to compare areal capacity and specific energy in a uniform, unambiguous manner. Gallagher et al.⁶ have discussed the shift of the current distribution due to salt depletion at the current collector. By scaling the transport equation, they derived an expression for the penetration depth of electrolyte transport
where ε_l is the porosity, τ_l is the tortuosity, D_l is the Li diffusion coefficient, t₊ is the transference number, c_l the salt concentration, and I_appl and F denote the current and the Faraday constant, respectively. The theoretical maximum areal capacity of an electrode can be expressed in terms of the specific capacity Q_s, the volume fraction ε_s, density ρ_s, and the electrode thickness L
with the definition of the penetration depth, the actual areal capacity can be expressed as
It should be noted that this description is only valid if the Lithium diffusion within the solid active material is not a limiting factor. The penetration depth and theoretical maximum areal capacity can be used to normalize the values obtained from the simulations.

Fig. 7 presents the normalized areal capacity and specific energy as a function of normalized thickness for all four electrode structures. The results show that the accessible capacity of electrodes without channels increases with the normalized thickness. In contrast, electrodes with equisized and narrow electrolyte channels maintain their accessible capacity up to a normalized thickness of approximately L/L_d = 12.5, while the capacity of electrodes with wide channels begins to decay at normalized thickness of approximately L/L_d = 7.5.

Fig. 7 Normalized areal capacity and specific energy of the four electrode structures as a function of normalized thickness at discharge rates of 1C, 2C, and 3C. The thickness was normalized with the electrolyte penetration depth introduced by Gallagher et al.¹¹ The curves for varying C rate collapse onto one curve for each electrode structure. The difference in electrochemical performance between hierarchical porous electrodes and electrodes without electrolyte channels can be clearly observed.

The benefit of the normalized presentation is even more apparent for the specific energy plots, where the curves collapse onto a single curve for each electrode structure. These plots clearly show the impact of the electrolyte channels on the specific energy that an electrode can deliver as a function of thickness. The performance of thick electrodes is substantially enhanced by the electrolyte channels compared with electrodes without channels. The specific energy of electrodes with wide channels decays slightly faster than the specific energy of electrodes with equisized and narrow channels, which show nearly the same characteristics.

Overall, our results show that the introduction of electrolyte channels thus can enhance the specific energy and reduce the fading at higher mass loading. However, this effect is limited by the size ratio of the electrolyte channels and the active material columns. If the electrolyte channels are too wide, the enhanced supply of Li salt through the macroporosity is offset by the increased thickness required to maintain mass loading, which leads to larger ohmic losses and reduces the accessible capacity.

5. Conclusions

We have performed numerical simulations of the discharge of hierarchical anisotropic electrode structures. The structures simulated can be experimentally produced by directional freeze casting followed by freeze drying and partial sintering. The resultant porous structure consists of continuous aligned parallel pore channels (macro-pores) through the thickness and the walls of the pore channels made up of isotropic micro-pores and the active electrode material. The important geometrical parameter is the ratio of the size of the macro pore channel (which gets filled with the electrolyte) and the wall thickness (w_elec : w_LTO). The model mimics freeze tape-cast electrodes, and the simulations were performed using experimentally measured materials properties of molybdenum-doped lithium titanate (Li₄Mo_xTi_5−xO₁₂) as active material and lithium hexafluorophosphate (LiPF₆) as electrolyte. We have assessed the rate capability at varying discharge rates and presented electrochemical performance characteristics in terms of capacity and specific energy. The numerical results confirm previous observations by Azami-Ghadkolai et al.³² that the introduction of anisotropic macroporosity can enhance the accessible capacity and improve the rate capability. Here, we investigated specifically the effect of the size ratio of the electrolyte channels and the active materials columns (w_elec : w_LTO). Our results indicate that the discharge characteristics are directly affected by the size of the electrolyte channels. While wide electrolyte channels limit the accessible capacity and specific energy due to the larger required thickness to maintain mass loading, electrodes with narrow channels show enhanced performance both in terms of accessible capacity and specific energy. The effect of the channel size has been rationalized in terms of the interplay of the salt concentration gradient and the current distribution in the active material during discharge. Electrodes with narrow channels discharge from the separator side and the current collector simultaneously such that only a small region in the center of the active material columns remains inaccessible.

While our results have been performed in two dimensions, the concentration profiles revealed that the main variation occurs in the direction of the current. It should be noted that our model geometries are a simplified representation of freeze tape-cast electrode structures and assume perfectly aligned and regular active materials columns. Real freeze tape-cast electrodes will deviate from perfect alignment and exhibit structural irregularities. While such irregularities would likely increase the tortuosity and reduce transport efficiency, this would affect any microstructure in a similar manner. Hence the relative effectiveness of electrolyte channels of varying width at mitigating salt depletion is expected to persist in practical freeze tape-cast electrodes. A future extension of our work is the use of digitized images of electrode structures to construct three-dimensional digital twins of real electrodes, which would allow us to probe the effect of structural irregularities and channel tortuosity on areal capacity and specific energy.

In conclusion, our results show that electrodes with hierarchical anisotropic porosity facilitate a higher mass loading than normal tape-cast electrodes and thus can maintain areal capacity and specific energy at higher discharge rates. The simulations further show that the size ratio of electrolyte channels and active material columns is a key factor in enhancing the electrochemical performance of electrodes. The simulations thus provide quantitative insight that can be leveraged in optimizing the design of hierarchical anisotropic porosity electrodes that provide high capacity and specific energy and enable rapid discharging and charging.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data for this article are available from Harvard Dataverse at https://doi.org/10.7910/DVN/A6SC20. The data analysis scripts of this article are available as an interactive capsule in the Code Ocean Open Science Library at https://doi.org/10.24433/CO.3737128.v1.

Acknowledgements

The authors would like to thank Dr. Milad Azami-Ghadkolai and Prof. Steve Creager for insightful discussions. We also thank Prof. Creager for generously providing access to the COMSOL Battery Design Module. Clemson University is acknowledged for a generous allotment of computer time on the Palmetto cluster. This work was supported by the National Science Foundation under award numbers OIA-1655740 and DMR-2119833, and the SC EPSCoR Program under award 20-GE05. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation nor the views of the SC EPSCoR Program.

References

1. U.S. Energy Information Administration, Monthly Energy Review, 2023.
2. G. E. Blomgren, The Development and Future of Lithium Ion Batteries, J. Electrochem. Soc., 2016, 164, A5019.
3. A. Manthiram, An Outlook on Lithium Ion Battery Technology, ACS Cent. Sci., 2017, 3, 1063.
4. J. T. Frith, M. J. Lacey and U. Ulissi, A non-academic perspective on the future of lithium-based batteries, Nat. Commun., 2023, 14, 420.
5. S. Kim, T. R. Tanim, E. J. Dufek, D. Scoffield, T. D. Pennington, K. L. Gering, A. M. Colclasure, W. Mai, A. Meintz and J. Bennett, Projecting Recent Advancements in Battery Technology to Next-Generation Electric Vehicles, Energy Technol., 2022, 10, 2200303.
6. K. G. Gallagher, et al., Optimizing Areal Capacities through Understanding the Limitations of Lithium-Ion Electrodes, J. Electrochem. Soc., 2015, 163, A138.
7. DE-FOA-0003162, Pioneering Railroad, Oceanic and Plane Electrification with 1 K Energy Storage Systems (Propel-1K), 2023.
8. T. Danner, M. Singh, S. Hein, J. Kaiser, H. Hahn and A. Latz, Thick electrodes for Li-ion batteries: a model based analysis, J. Power Sources, 2016, 334, 191.
9. B. Scrosati, J. Garche and W. Tillmetz, Advances in Battery Technologies for Electric Vehicles, Woodhead Publishing, 2015.
10. R. Zhao, J. Liu and J. Gu, The effects of electrode thickness on the electrochemical and thermal characteristics of lithium ion battery, Appl. Energy, 2015, 139, 220.
11. C. Fongy, A.-C. Gaillot, S. Jouanneau, D. Guyomard and B. Lestriez, Ionic vs Electronic Power Limitations and Analysis of the Fraction of Wired Grains in LiFePO₄ Composite Electrodes, J. Electrochem. Soc., 2010, 157, A885.
12. A. Smith, P. Stüble, L. Leuthner, A. Hofmann, F. Jeschull and L. Mereacre, Potential and Limitations of Research Battery Cell Types for Electrochemical Data Acquisition, Batteries Supercaps, 2023, 6, e202300080.
13. Y. Kim, M. Kim, T. Lee, E. Kim, M. An, J. Park, J. Cho and Y. Son, Investigation of mass loading of cathode materials for high energy lithium-ion batteries, Electrochem. Commun., 2023, 147, 107437.
14. A. M. Boyce, X. Lu, D. J. L. Brett and P. R. Shearing, Exploring the influence of porosity and thickness on lithium-ion battery electrodes using an image-based model, J. Power Sources, 2022, 542, 231779.
15. T. Beuse, M. Fingerle, C. Wagner, M. Winter and M. Börner, Comprehensive Insights into the Porosity of Lithium-Ion Battery Electrodes: A Comparative Study on Positive Electrodes Based on LiNi_0.6Mn_0.2Co_0.2O₂ (NMC622), Batteries, 2021, 7, 70.
16. J. Landesfeind, A. Eldiven and H. A. Gasteiger, Influence of the Binder on Lithium Ion Battery Electrode Tortuosity and Performance, J. Electrochem. Soc., 2018, 165, A1122.
17. Y. Dai and V. Srinivasan, On Graded Electrode Porosity as a Design Tool for Improving the Energy Density of Batteries, J. Electrochem. Soc., 2015, 163, A406.
18. M. A. Ghadkolai, S. Creager, J. Nanda and R. K. Bordia, Freeze Tape Cast Thick Mo Doped Li₄Ti₅O₁₂ Electrodes for Lithium-Ion Batteries, J. Electrochem. Soc., 2017, 164, A2603.
19. B. Delattre, R. Amin, J. Sander, J. D. Coninck, A. P. Tomsia and Y.-M. Chiang, Impact of Pore Tortuosity on Electrode Kinetics in Lithium Battery Electrodes: Study in Directionally Freeze-Cast LiNi_0.8Co_0.15Al_0.05O₂ (NCA), J. Electrochem. Soc., 2018, 165, A388.
20. S. Behr, R. Amin, Y.-M. Chiang and A. Tomsia, Highly-Structured, Additive-Free Lithium-Ion Cathodes by Freeze-Casting Technology, Ceram. Forum Int., 2015, 92, 39.
21. J. Newman and W. Tiedemann, Porous-electrode theory with battery applications, AIChE J., 1975, 21, 25.
22. T. F. Fuller, M. Doyle and J. Newman, Simulation and Optimization of the Dual Lithium Ion Insertion Cell, J. Electrochem. Soc., 1994, 141, 1.
23. K. E. Thomas, J. Newman and R. M. Darling, Mathematical Modeling of Lithium Batteries, in Advances in Lithium-Ion Batteries, ed. W. A. van Schalkwijk and B. Scrosati, Springer US, Boston, MA, 2002, pp. 345–392.
24. J. Newman and K. E. Thomas-Alyea, Electrochemical Systems, John Wiley & Sons, 2012.
25. P. D. Vidts and R. E. White, Governing Equations for Transport in Porous Electrodes, J. Electrochem. Soc., 1997, 144, 1343.
26. G. G. Botte, V. R. Subramanian and R. E. White, Mathematical modeling of secondary lithium batteries, Electrochim. Acta, 2000, 45, 2595.
27. E. García, Y.-M. Chiang, W. Carter, P. Limthongkul and C. Bishop, Microstructural Modeling and Design of Rechargeable Lithium-Ion Batteries, J. Electrochem. Soc., 2005, 152, A255.
28. D. Danilov and P. H. L. Notten, Mathematical modelling of ionic transport in the electrolyte of Li-ion batteries, Electrochim. Acta, 2008, 53, 5569.
29. A. Latz and J. Zausch, Thermodynamic consistent transport theory of Li-ion batteries, J. Power Sources, 2011, 196, 3296.
30. A. Latz and J. Zausch, Multiscale modeling of lithium ion batteries: thermal aspects, Beilstein J. Nanotechnol., 2015, 6, 987.
31. S. Allu, S. Kalnaus, S. Simunovic, J. Nanda, J. A. Turner and S. Pannala, A three-dimensional meso-macroscopic model for Li-Ion intercalation batteries, J. Power Sources, 2016, 325, 42.
32. M. Azami-Ghadkolai, M. Yousefi, S. Allu, S. Creager and R. Bordia, Effect of isotropic and anisotropic porous microstructure on electrochemical performance of Li ion battery cathodes: an experimental and computational study, J. Power Sources, 2020, 474, 228490.
33. R. B. Smith and M. Z. Bazant, Multiphase Porous Electrode Theory, J. Electrochem. Soc., 2017, 164, E3291.
34. A. G. Kashkooli, S. Farhad, D. U. Lee, K. Feng, S. Litster, S. K. Babu, L. Zhu and Z. Chen, Multiscale modeling of lithium-ion battery electrodes based on nano-scale X-ray computed tomography, J. Power Sources, 2016, 307, 496.
35. R. de Levie, On porous electrodes in electrolyte solutions: I. Capacitance effects, Electrochim. Acta, 1963, 8, 751.
36. A. G. Kashkooli, G. Lui, S. Farhad, D. U. Lee, K. Feng, A. Yu and Z. Chen, Nano-particle size effect on the performance of Li₄Ti₅O₁₂ spinel, Electrochim. Acta, 2016, 196, 33.
37. M. Rashid, A. Sahoo, A. Gupta and Y. Sharma, Numerical modelling of transport limitations in lithium titanate anodes, Electrochim. Acta, 2018, 283, 313.
38. C. L. Cobb and M. Blanco, Modeling mass and density distribution effects on the performance of co-extruded electrodes for high energy density lithium-ion batteries, J. Power Sources, 2014, 249, 357.
39. M. Doyle, T. F. Fuller and J. Newman, Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell, J. Electrochem. Soc., 1993, 140, 1526.
40. M. Doyle, J. Newman, A. S. Gozdz, C. N. Schmutz and J.-M. Tarascon, Comparison of Modeling Predictions with Experimental Data from Plastic Lithium Ion Cells, J. Electrochem. Soc., 1996, 143, 1890.
41. W. Lai and F. Ciucci, Mathematical modeling of porous battery electrodes—Revisit of Newman's model, Electrochim. Acta, 2011, 56, 4369.
42. J. Huang and Y. Zhang, Essays on Conceptual Electrochemistry: I. Bridging Open-Circuit Voltage of Electrochemical Cells and Charge Distribution at Electrode–Electrolyte Interfaces, Front. Chem., 2022, 10, 938064.
43. J. Ni and C. Beckermann, A volume-averaged two-phase model for transport phenomena during solidification, Metall. Trans. B, 1991, 22, 349.

---