Li+ concentration waves in a liquid electrolyte of Li-ion batteries with porous graphite-based electrodes

Electrolyte solutions function as ionic conductors in Li-ion batteries and inevitably induce concentration gradients during battery operation. It is shown that in addition to these concentration gradients, very specific Li + concentration waves in the electrolyte are formed in graphite-based porous electrode/Li cells. This phenomenon has been investigated by both simulations and experiments. From the simulations, it has been concluded that the occurring Li + concentration waves in the electrolyte vary with position and time. Such waves originate from the fluctuations of the reaction distribution inside the porous electrode and depend on both the thermodynamics (open-circuit voltage, OCV) and kinetics (charge transfer reaction heterogeneity). Li + concentration waves occurring inside the separator region are directly related to the battery output voltage at low current applications. A four-electrode device is used to validate the electrolyte concentration waves experimentally. The electric potential differences between the reference electrodes and counter electrode show regular fluctuations, demonstrating the existence of concentration waves in the electrolyte. The simultaneous appearance of the fluctuations in the potential differences and the transitions from plateaus to slopes in the battery output voltage illustrates the dependency of Li + concentration waves on the thermodynamics and kinetics of the electrochemical reactions.


Introduction
Since their successful commercialization in 1991, Li-ion batteries (LIBs) have been extensively applied as energy storage devices for portable electronics and (hybrid) electric vehicles due to the highly favorable combination of reduced cost and impressive battery performance [1][2][3] . Typical LIBs consist of two porous electrodes and a porous separator impregnated by an electrolyte. The most adopted electrolytes for LIBs are non-aqueous solutions, in which lithium hexafluorophosphate (LiPF 6 ) salt is dissolved in the mixture of solvents. These solvents include organic esters and ethers [ 4 , 5 ], i.e., ethylene carbonate (EC), dimethyl carbonate (DMC), diethyl carbonate (DEC), dimethoxyethane (DME), etc . LiPF 6 in the mixture of solvents dissociates into solvated Li + and PF 6 − [ 6 , 7 ]. The electrochemical charge-transfer reactions occur at the interfaces between electrodes and electrolytes during operation. These reactions are accompanied by the movement of electrons and ions [8] . Electrons are transported through the external circuit to/from the electrode's active material. Ions move through the porous electrodes and the electrolyte and are driven by diffusion and migration. Due to ionic transport limitations, concentration gradients in the electrolyte are de-veloped [9][10][11] . These gradients are important sources of polarization [ 11 , 12 ]. It is commonly accepted that the charge-transfer reaction also depends on the salt concentration in the electrolyte [13][14][15] . Therefore, understanding the electrolyte behavior will help unravel the essence of the electrochemical reactions and may shed further light on the optimal design and operation of LIBs.
Simulations and experiments have been extensively applied to investigate the electrolyte concentration during battery operation. Physicsbased electrochemical models are frequently adopted to perform these simulations [16][17][18][19][20][21][22][23][24] . An important output of such models is the Li + concentration distribution in the electrolyte across the battery thickness as a function of time and position [ 10 , 16-18 , 25 ]. According to simulations [ 12 , 26-28 ] and experimental investigations [ 9 , 29 ], the Li + concentration in the electrolyte at the delithiated (porous) electrode is larger than that at the lithiated (porous) electrode. This phenomenon is explained by the movements of Li-ions across the electrolyte/electrode interfaces. During lithiation of the electrode material, Li-ions are extracted from the electrolyte, lowering the salt concentration in the solution. An inverse behavior is observed during delithiation. Ionic transport is quantified by transport parameters, such as the diffusion coefficient and transference number of cation in the electrolyte [30] . Smaller diffusion coefficients and/or lower cation transference numbers lead to larger electrolyte concentration gradients. Moreover, these two parameters are mathematically related to each other and are influenced by concentration and temperature [ 30 , 31 ]. In addition, other factors will also influence the electrolyte concentration distribution, such as the applied current density, porosities, thickness of the electrode and separator, electrolyte convection, etc .
Several experimental techniques have been developed to in situ determine the electrolyte concentration distribution during battery operation, such as X-ray based imaging [ 9 , 32 ], nuclear magnetic resonance (NMR), magnetic resonance imaging (MRI) [ 29 , 33 , 34 ], Raman microscopy [35][36][37] , and neutron depth profiling (NDP) [38][39][40] . Although these methods provide valuable information and all show the distribution of Li-ions inside the electrolyte, these methods are all based on highly sophisticated equipment. High noise levels sometimes may cover interesting but subtle signals. Electrochemical methods are a good alternative due to their easy accessibility, high sensitivity, and real-time response. Electrical potential differences measured between the (reference) electrodes are related to the electrochemical potential differences of the electrodes in the solution and depend on the salt concentration in the electrolyte [30] . According to this relation, the variations of the Li-ion concentration in the electrolyte can be extracted [41][42][43] .
In the present work, batteries composed of graphite-based porous electrodes and Li metallic counter electrodes (C/Li) have been selected as a research subject to reduce the complexity of complete cells consisting of two porous electrodes. The Li-ion concentration in the electrolyte during galvanostatic operation is systematically analyzed from both simulation and experimental perspectives. Traditionally, it is believed that the Li + concentration gradients in the electrolyte are stabilized after entering the steady state of the charge transfer reaction. Our results reveal that the Li + concentration does not completely stabilize but shows remarkable waves during the operation of C/Li batteries. These waves are created by fluctuations of the reaction distribution in the porous graphite-based electrode, which strongly depend on the reaction thermodynamics and kinetics.

Experimental
Coin-type cells: Commercial 18650-type cylindrical batteries manufactured by Tianjin Lishen Battery Co., Ltd were dismantled in an argonfilled glove box, and pieces of the anodes were taken out. The active material in the anode is graphite with a small amount of silicon [44] . It is therefore called 'graphite-based electrode' in this manuscript. Before the electrochemical measurements, the active material on one side of the double-coated graphite-based electrodes was carefully removed with the help of acetone and a sharp blade. Then, the remainder was cleaned with dust-free tissues immersed in acetone. This procedure ensures that the active material is completely removed and scratching the current collector is only minor. Then the intact parts of graphite-based electrodes were selected and cut into discs with a diameter of 14 mm. 2032-type coin cells were subsequently assembled using the as-prepared graphite-based electrodes as working electrodes and Li metal foil as counter electrodes. A 2400-type Celgard separator (25 μm thick) and 1 M LiPF 6 electrolyte in a solvent mixture of EC:DMC:DEC with a 1:1:1 volume ratio were used.
The electrochemical properties of the assembled cells were measured by an automated battery cycler (Neware) in the voltage range of 0.01-2 V at 25°C. Before the test, all cells were equilibrated for 12 h and then activated for four cycles at 0.2 C-rate (1C = 7 mA) in the constantcurrent charging mode (delithiation of the graphite-based electrode) and constant-current constant-voltage (CCCV) discharging mode (lithiation of graphite-based electrode). The voltage and cut-off current in the CVmode were 0.01 V and 0.04C, respectively. Characterization cycles were then performed in the CCCV-mode, using a constant discharge current (0.2C), followed by a resting period of 30 mins and a set of constant charging currents (0.04, 0.1, 0.2, 0.5, 1.0, and 1.4C) applied in sequence in the subsequent cycles. The 30-minute resting period was selected because the voltage showed the maximum changes within this period for discharge. The voltage showed only a minor increase after a longer resting time. In the simulations, the potential curve obtained at 0.01 C-rate and 25°C was used as the (pseudo) open-circuit voltage (OCV) curve. Other researchers frequently adopted the potential curve at low C-rates as OCV in the simulations and experiments [45][46][47] .
Four-electrode cells: For the assembly of the reference electrode cells, the as-prepared graphite-based electrodes were cut into discs with a diameter of 18 mm and used as the working electrode (WE) in the fourelectrode EL-Cell (PAT-Cell-TwinRef). Li metal with a diameter of 18 mm was used as the counter electrode (CE). Separators FP-5S were used, with 21.6 mm in diameter and 220 μm in thickness. Two ring-shaped Li-reference electrodes, denoted as RE1 and RE2, were assembled at the two sides of the separator. 100 μL of the electrolyte mentioned above was used.
The electrochemical properties of the four-electrode EL-Cell were measured with a VMP 300 potentiostat (Bio-Logic) in the voltage range of 0.005-1.6 V under temperature control at 25°C. Two channels were used to record the potential differences among the four electrodes. One channel provided the power and recorded the mutual potential differences among WE, CE, and RE2. Another channel recorded the potential difference between RE1 and RE2. Before the test, the as-prepared cells were put to rest for 20 h and then activated for three cycles at 0.05 C-rate (1C = 12 mA). Subsequently, 0.05, 0.1, 0.2, and 0.4C charging currents were applied in sequence, followed by a discharging current of 0.05C. The cells were also subjected to 0.05, 0.1, and 0.2 C-rate for discharging and 0.05 C-rate for charging. Between the charging and discharging, a resting period of 6 h was adopted.
For assembling the four-electrode EL-Cell with LiFePO 4 (LFP) as WE, LFP electrodes were obtained by dismantling a commercial 26650-type LFP/graphite cylindrical battery and cutting the LFP electrode into discs with a diameter of 18 mm. The other parts were the same as used in the graphite-based four-electrode EL-Cell. For the measurements of the LFP-based four-electrode EL-Cell, the regime was kept similar to the cell with graphite-based WE. The voltage range was set at 2.8-4.2 V, and the current took values 0.1, 0.2, 0.4, and 0.8 C (1C = 4 mA) to keep the current density similar to that of the graphite-based EL-Cell. Note that the 1 C-rate for coin-type cells with graphite-based WE equals 7 mA, which has been confirmed by multiple cell tests. The corresponding 1Crate for the EL-Cell with graphite-based WE is 12 mA, and 4 mA for LFP WE's.
Simulations: Simulations were performed using the Matlab R2018b software package. The finite volume method [23] was used to discretize the partial differential equations (PDE) listed in Table S1 into systems of algebraic equations. A numerically efficient method [48] was adopted to solve the full-order pseudo-two-dimensional (P2D) model, in which the reaction distribution was first optimized and subsequently used as input for the calculations of the Li concentration in the electrolyte and solid. Parameter definitions and parameter values are listed in Table S2 and S3, respectively. The temperature in the simulations was set to 25°C. It is the same as the temperature in all experiments.   Table S3. A pseudo-two-dimensional (P2D) model developed by Newman et al. [ 16 , 17 ] describes the physical and electrochemical processes inside this battery. The governing equations and boundary conditions are listed in Table S1. The parameter definitions and parameter values are given in Table S2 and S3, respectively. Fig. 1 b shows the comparison between the simulated (black line) and experimental (symbols) voltage curves of a coin-type cell at 0.1C charging current (1C = 7 mA). Note that the charging process denotes the delithiation of the graphite-based electrode. Good agreement can be found between the simulation and experiment. A comparison of the voltage curves at other C-rates can be found elsewhere [44] . The red curve in Fig. 1 b gives the open-circuit voltage (OCV) curve. The plateaus and sloping regions in the OCV curve of the graphite electrode are related to the two-phase (stage) coexistence regions and solid-solution region [49][50][51] , respectively. Fig. 1 c shows the simulated profiles of Li + concentration in the electrolyte ( 2 ) as a function of normalized position and time. Typical 2 distributions are observed during delithiation, implying high 2 inside the porous electrode and low 2 concentrations in the separator region [ 17 , 21 ]. Strikingly, 2 shows clear fluctuated waves in both the porous electrode and separator regions. Such fluctuations of 2 have been observed before [ 19 , 23 , 52 ]. However, rather little attention is paid to this striking phenomenon.

Concentration waves in the electrolyte
From the governing equations, only mass transport is responsible for the evolution of 2 in the separator region (Equation S3 in Table  S1) while both mass transport and the reaction distribution ( C ) influence 2 in the porous electrode region (Equation S9 in Table S1). C is calculated from the Butler-Volmer equation and can quantify the distribution of charge-transfer reaction rate inside the porous electrode. In other publications, C can also be called the pore wall flux [ 16 , 17 ] or reaction flux distribution [53][54][55] .
The parameters for mass transport, i.e., the diffusion coefficient of Li + in the electrolyte ( 2 ) and the transference number of Li + in the electrolyte ( + ), are kept constant in the present simulations. Intuitively, C causes the waves in 2 . Fig. 1 d shows C inside the porous graphitebased electrode at 0.1C charging as a function of the normalized position and time together with the projected two-dimension (2D) contour image on the top. It can indeed be seen that C is not uniform during charging. At very short timescales, close to 0 h, C shows a larger distribution near the SC interface. Between 0 and 8 h, three main peaks in C initiate at the SC interface. These peaks propagate progressively from the SC to the CC interface, forming waves. From 8 h to the end of charging, a relatively uniform C can be found. The three main waves in C agree well with the three main waves observed in 2 ( Fig. 1 c). The relatively uniform C after 8 h also corresponds to the minor fluctuations in 2 .
To confirm the dependence of 2 on C , Fig. S1a shows 2 with a uniform reaction distribution ( C ) across the porous electrode ( Fig. S1b), given by Such an assumption implies that the reaction is uniformly distributed across the porous electrode, leading to a constant C . The derivation of this equation can be found in Section II of the Supporting Information.
A similar equation can also be found in the literature [ 44 , 56 , 57 ]. The model based on this simplification is called average model. In this case

Reaction distribution waves
Reaction distribution C is influenced by the exchange current density ( 0 C ), the charge-transfer overpotential ( C ) and the charge-transfer coefficient ( ) of the graphite-based electrode according to Equation S6 in Table S1. In the present work, is taken as a constant equal to 0.5. Consequently, only 0 C and C influence the reaction distribution C . 0 C is a function of the Li + concentration at the surface of solid particles ( 1 ) and the Li + concentration in the electrolyte ( 2 ) as expressed by Equation S8. C is a function of the electrical potential in solid ( Φ 1 ), the electrical potential in the electrolyte ( Φ 2 ), and the surface potential of the solid particles ( C ) as described by Equation S7. Note that C can be interpreted as the electrode OCV ( C ) evaluated with Li concentration at the particle surface 1 , i.e. C = C ( 1 ) . To clarify the relations, Fig. 2 gives the profiles of these variables with respect to the normalized position at three typical moments of time, namely 0.13, 1.33, and 3.33 h. These three moments are indicated by the vertical arrows and the (partly hidden) white lines in Fig. 1 d. The capacities corresponding to these moments are written near the arrows in Fig. 1 b. These moments represent the time before the wave formation, the wave initiation, and the end of the wave. Fig. 2 a-c shows the development of C at these three moments. In Fig. 2 a, C has the highest value near the SC interface. In Fig. 2 b, a mound-shaped C is formed with the maximum near the SC interface at 1.33 h. This maximum slowly propagates towards the CC interface at 3.33 h ( Fig. 2 c).
shows the corresponding development of 0 C and C inside the porous electrode at these moments. The changing trends of C in Fig. 2 d-f are similar to those of C in Fig. 2 a-c, indicating that they are strongly correlated. The peaks of C locate at the same position with those of C , as marked by the blue shade areas. 0 C shows an attenuated peak in Fig. 2 f, but this peak does not coincide with that of C . In the present paper, 0 C is a function of 1 and 2 as expressed by Equation S8. Fig. S2 shows 0 C with respect to 1 keeping 2 equal to 1000 mol/m 3 . One can see that 0 C has a mound shape. The attenuated peak of 0 C in Fig. 2 f is, therefore, due to the influence of changing 1 . To eliminate the influence of the changing 1 and 2 , a constant 0 C was used for the simulations of Fig. S3 and S4. In these figures, waves still can be seen in C and C . Therefore, the concentration-dependent 0 C does not result in the multiple peaks in C . The varied charge-transfer reaction overpotential C causes the multiple peaks in C .
To explore the relationship between C and C , Fig. S5 shows the various components contributing to C (Equation S7), namely Φ 1 , Φ 2 , and C at the three indicated moments. To help in understanding, Fig. 2 g-i shows Φ 1 − Φ 2 and C . Note that the difference between Φ 1 − Φ 2 and C equals to C (red-shaded areas). The values of Φ 1 − Φ 2 in Fig. 2 g-i decrease smoothly from the SC to the CC interface due to a large effective electronic conductivity ( C ) of the solid and a small effective ionic conductivity ( C ) of the solution in the porous electrode region, and the Li + concentration gradient in the electrolyte. Φ 1 − Φ 2 increases continuously in time due to Li extraction from the material (for the separate de-velopment of Φ 1 and Φ 2 , please refer to Fig. S5). C in Fig. 2 g-i shows an interesting distribution along the normalized position. A relatively flat profile is observed in Fig. 2 g. A bending area of C initiates near the SC interface in Fig. 2 h, marked by the blue shaded area. In Fig. 2 i, the bending area appears near the CC interface (blue shaded area). These bending areas of C with respect to the position ( Fig. 2 h-i) cause the peaks of C ( Fig. 2 e-f) with the help of Φ 1 and Φ 2 , which finally result in the peaks of C ( Fig. 2 b-c), as shown by the blue shaded areas.
To investigate the bending area of C along with the position, Movie 1 in the Supporting Information illustrates the dynamic changes of the variables mentioned above at the full duration of the galvanostatic charging. The following quantities are plotted. The battery voltage and OCV are plotted vs. the transferred charge in subplot (a). C as a function of 1 and the corresponding derivative are in subplot (b). Subplot (c) contains C as a function of the position. Subplot (d) shows overpotential C . Finally, subplot (e) gives C with respect to position, and (f) contains the selected C at the SC and CC interfaces accordingly. All these quantities change with charging time. Note that C in (b) is plotted with respect to the surface Li concentration of particles 1 , which is a function of time and normalized position. Consequently, C can also be plotted vs. time and normalized position. These two plots are given in Fig. S6 and (c) in Movie 1, respectively. The blue, black, and red dots in (b-f) represent the particle at the SC interface, in the middle of the porous electrode, and at the CC interface. At the very beginning (0 h), C is determined only by the battery parameters, i.e. C , C , 0 C , , and [ 55 , 58 ], etc. , without the influence of C . For the given set of parameters C is dominant near the SC interface, causing 1 at this region decreases faster than other regions. The decreased 1 tends to influence the distribution of C but only throughout the OCV curve. At the plateau of C vs. 1 (between 0 and 1.1 h in (b)), C shows a larger value near the SC interface (e) due to a relatively flat distribution of C vs. positon in (c).
Because 0 C is larger near the SC interface (see Fig. 2 d-e), C still shows a dominant distribution near the SC interface (e). That, in turn, causes 1 near the SC interface decreases even faster than that in other regions (b).
Subsequently, the particle at the SC interface (blue dot in (b)) first enters the sloping region of C (vs. 1 ), resulting in a bending area of C vs. position near the SC interface (c). The bending area of C vs. position causes the disturbances of C , forming peak (local maximum) at the SC interface (d), which further causes the peak formation of C at the SC interface (e-f). After that, various particles located between the SC and CC interface (between the blue and red dot in (b)) subsequently enter the sloping region from the plateau of C (vs. 1 ), which leads to the movement of the bending area of C (vs. position) in (c). This movement of the bending area of C results in the peak propagation of C and C from the SC to the CC interface in (d-f), finally forming the first wave as observed in Fig. 1 d. The second and third waves of C in Fig. 1 d are all related to the movement of the bending area of C vs. position.
From Fig. 2 and Movie 1, several conclusions can be drawn. C tends to become uniform across the normalized position of the porous electrode at long-time scales. The plateaus of C (vs. 1 ) make C enter the uniform state slowly, causing a wide distribution of 1 and further leading to large inhomogeneous utilization of the porous electrode. The sloping regions of C (vs. 1 ) make C enter the uniform state quickly, causing a relatively narrow distribution of 1 and further leading to slightly inhomogeneous utilization of the porous electrode. The transitions from the plateaus to the sloping regions of C (vs. 1 ) generate the waves of C . These waves mitigate the large inhomogeneous utilization brought by the plateaus of C .
Specifically, two factors influence the C wave behavior. The first one is transitions from plateaus to slopes (steps) in C curve, which cause steps in C . C curve is thermodynamically determined by the material structure and the nature of reaction [59][60][61] . The second factor is the nonuniform distribution of 1 along with the normalized position, which can also be interpreted as the nonuniform electrode utilization. The electrochemical reaction heterogeneity plays a vital role in determining the nonuniform electrode utilization. These two factors are further denoted as the thermodynamic and kinetic factors and are discussed below in detail.
C is the battery OCV curve (or equilibrium potential), which is thermodynamically determined by the electrochemical potential difference between the cathode and anode [59][60][61] . The properties of C directly influence surface potential of solid particles C as a function of 1 . From Fig. 2 and Movie 1, transitions from plateaus to slopes (steps) in C cause the peak formation and propagation in C . Three major steps in C , heritated from C , result in three major waves of C , which can be observed in Fig. 1 d. To explore the influence of steps in the OCV curve on C and 2 , the experimental C , shown in Fig. 1 b, is numerically modified into curves with two (case 1), one (case 2), and none (case 3) steps, respectively, as shown in Fig. 3 a-c. The removed steps are replaced with a smoothly increasing curve. The experimental OCV curve is also shown in Fig. 3 a-c for comparison (grey curves). The simulated C and 2 for these three modified OCV-curve cases are shown in Fig. 3 d-f and Fig. 3 g-i, respectively. In case 1 ( Fig. 3 a), the step at the SoC range of 0.7-0.8 is smoothed, and the steps at the SoC of 0.2-0.4 and 0.4-0.6 remain. The corresponding C ( Fig. 3 d) and 2 ( Fig. 3 g) show two waves. The third wave observed in Fig. 1 c-d disappears. In case 2 ( Fig. 3 b), the steps at the SoC range of 0.7-0.8 and 0.4-0.6 are smoothed, while the step at SOC of 0.2-0.4 remains. Correspondingly, only one wave can be observed in C ( Fig. 3 e) and 2 ( Fig. 3 f), and the other two waves observed in Fig. 1 c-d disappear. In case 3 ( Fig. 3 c), the OCV curve is completely smoothed without any steps. C ( Fig. 3 f) and 2 ( Fig. 3 i) consequently show smooth profiles without any waves.
Another essential factor causing the waves of C is the nonuniform electrode utilization, resulting from the heterogeneity of electrochemical reaction. Essentially, electrode utilization refers to the particle delithiation state as a function of the normalized position. The particle delithiation state is affected by the amount of lithium extracted, which is determined by the integral of C with respect to time. In other words, the nonuniform C leads to nonuniform electrode utilization. Several parameters influence the distribution of C , i.e. C , C , 0 C , , , etc. , at the short-time scales [ 55 , 58 ], and gradient in 1 and 2 and the OCV curve at the long-time scales. All these parameters and the related factors will affect the heterogeneity of electrochemical reaction, and further influence the waves of C . Fig. 4 a and d show the results of the simulations for C and 2 with enlarged effective ionic conductivity C and electrolyte diffusion coefficient C 2 . An experimental OCV curve was used for these figures. The larger C together with a large effective electronic conductivity C (beyond the saturation value, which is determined by the electrode parameters) makes C uniform at short-time scales [58] . Enlarged C 2 eliminates the influence of electrolyte concentration gradient on Φ 2 at longer time scales (Equation S10). Such a combination of parameters makes C uniform at both the short-and long-time scales, as shown in Fig. 4 a. In this case, C shows a smooth and flat distribution during the entire charging process without any waves. The minor variations are due to minor electrolyte concentration gradients. Fig. 4 b and e simulate the case of using a small C combined with a large C . For such combination of parameters [58] , C becomes dominant near the CC interface before the wave generation ( Fig. 4 b). In this case, the waves in C can also be found, but these waves now propagate from the CC to the SC interface, in the opposite direction of Fig. 1 d. Fig. 4 c and f show the case using a small C and a small C . In this case C is dominant at both the SC and CC interfaces at short-time scales ( Fig. 4 c) [58] . For this specific combination of parameters, the waves of C originate from both the SC and CC interfaces and propagate to the middle of the porous electrode, where they converge.
In summary, both the OCV curve (thermodynamic factor) and the electrochemical reaction heterogeneity (kinetic factor) are essential for the generation and propagation of waves in C . The reaction heterogeneity causes the nonuniform electrode utilization across the battery. This nonuniform electrode utilization tends to be mitigated when the electrode reaction is going through transitions from the plateaus to the sloping regions in the OCV curve, forming waves of C .  Fig. 5 a-b. The peak of the Li + concentration in the electrolyte ( 2 ) inside the porous electrode region clearly propagates from the SC to the CC interface as indicated by red arrows in Fig. 5 b. Such behavior is in line with the development of the three C waves in Fig. 5 a. Low concentrations in 2 emerge in the separator region when 2 is high at the CC interface (blue-shaded areas in Fig. 5 b). This behavior is attributed to the mass conservation of 2 across the battery, i.e. ∫ 0 2 = 0 2 , where 0 2 denotes the initial Li + concentration in the electrolyte and the total pore volume inside the battery, including the separator and the graphite-based electrode. As a consequence, when high concentrations in 2 appear as peaks at one interface, low concentrations in 2 inevitably appear as valleys near the other interface. The waves of C inside the porous electrode are responsible for the 2 waves in both the porous electrode and separator regions. When no waves exist in C ( Fig. 3 f and S1b), 2 shows smooth profiles in both the porous electrode and separator region ( Fig. 3 i and S1a).

Wave-dependence of electrolyte concentration on reaction distribution
Ionic transport in the electrolyte is also essential for 2 behavior. Figs. S7 and S8 show C and 2 with various electrolyte diffusion coefficient ( 2 ) and transference number of Li + in the electrolyte ( + ). In these figures, C all shows three major waves propagating from the SC to the CC interface. Minor changes, i.e., wave position, wave width, and height,   Table S3. The effective electrolyte diffusion coefficient ( C 2 ), effective ionic ( C ) and electronic conductivity ( C ) in the porous graphite electrode region should be modified with respect to the porosity and Bruggeman coefficient, according to  can be observed compared to Fig. 5 a due to the influences of 2 . However, 2 shows considerable changes with different 2 or + . Note that 2 -gradients refer to the Li + concentration distribution along with the position, as indicated by the black vertical lines in Fig. S7-S8. 2 -waves refer to the Li + concentration fluctuations with respect to time and position, as indicated by the red arrows in Fig. 5 b. Large values of 2 or + effectively suppress the 2 -gradients ( Fig. S7b and Fig. S8b). Small values of 2 or + lead to larger gradients of 2 (Figs. S7d and S8d). However, 2 -waves still exist with different 2 or + even though they become less prominent with larger 2 or + . That is because 2 -waves are generated by the waves of C , which is minorly influenced by the ionic transport parameters. Therefore, the evolution of 2 in the porous electrode and separator regions is controlled by both C inside the porous electrode and the ionic transport in the electrolyte. The waves in 2 originate from the waves of C inside the porous electrode. The gradients in 2 are attributed to the ionic transport limitation in the electrolyte.
Another factor that influences the 2 -gradients and waves is the Crate. Fig. S9 shows the voltage comparison and Li + concentration waves at various C-rates using the parameters listed in Table S3. Note that the x -axis has been replaced by transferred charge to make them comparable. At low C-rate, the Li + concentration waves are prominent. At a high C-rate (larger than 1 C), Li + concentration waves can still be observed, but they are not as obvious as at low C-rates. The reason for such a behavior is that large Li + concentration gradients across the thickness of electrode material develop during high C-rates operation. Parts of the electrodes at different depths create concentration waves at different moments in time, smoothing the total concentration wave.
As discussed above, the transitions from the flat to the sloping parts (steps) in C (vs. 1 ) coincide with the waves in C . However, C is not a measurable quantity, making it difficult to validate the simulations. To relate C to easily measurable quantities, for example, the battery voltage ( ), an expression of (Equation S15) has been derived in Section III of the Supporting Information based on equations in Table  S1. It can be seen that is expressed as the summation of the surface potential of electrode particles at the CC interface ( C | = ) and three additional overpotential terms. Note that the overpotential terms are very small at low C-rate operation. Figs. 5 c and S10 present , C | = and battery OCV curve at 0.1 and 0.04C, accordingly. It can be seen that (black curve) behaves similarly to C | = (red curve) since the overpotentials are small for such low C-rates. Even though the fluctuating overpotential components [44] let the plateaus and slopes in slightly shift with respect to those found for C | = , still reflects the general features of C | = . As depicted in Fig. 5 b, the peaks of 2 at the CC interface correspond to the transitions from the flat to the sloping regions (steps) in C | = . The valleys of 2 at the separator region nicely align with the 2 -peaks at the CC interface (see blue-shaded areas in Fig. 5 b-c) due to mass conservation in the electrolyte. It can therefore be concluded that the valleys of 2 at the separator region also appear near the steps of , as shown by the blue-shaded areas in Fig. 5 b-c.

Experimental validation
The P2D-model adopted in this work combines a macroscale dimension (along with the battery thickness) and a microscale dimension (along with the particle radius). The waves of 2 are generated in the macroscale dimension. It is not easy to measure 2 inside porous electrodes. However, 2 across the separator can conveniently be indicated by the potential difference measured with reference electrodes [ 41 , 62 ]. The waves of 2 at the separator region behave oppositely with those at the CC interface and are related to , as explained in Fig. 5 . Therefore, 2 across the separator is selected to verify the Li + concentration waves in the electrolyte. An EL-Cell with two Li ring reference electrodes (PAT-Cell-TwinRef, Germany) is adopted. Fig. 6 a shows the schematic of this cell. A graphite-based electrode used in Fig. 1 -5 is assembled as the working electrode (WE). A metallic Li electrode is the counter electrode (CE). Two stainless steel rings coated with metallic Li are used as the reference electrodes RE1 and RE2, assembled at two sides of an electrolyteimmersed separator. Note that RE1 locates at the WE side. The potential differences among CE, RE1, and RE2 are mathematically derived in Equation S16-S28 in Section IV of the Supporting Information. Apparently, the potential differences between CE and RE2 (E CE -E RE2 ) and between RE1 and RE2 (E RE1 -E RE2 ) are all a function of 2 across the separator region. Fig. 6 b shows the potential differences between the various electrodes of the EL-Cell during charging and discharging at 0.05C (0.24 mA/cm 2 ) followed by 6 h resting periods. The top panel shows the potential difference (light blue curve) between the WE and CE (E WE -E CE ) and the derivative of this potential difference dV/dt (light red curve). Typical plateaus and slopes can be discerned in the voltage curve of E WE -E CE . The positions of the transitions from the plateaus to slopes (steps) are emphasized by the peaks in the potential derivative curve. The intermediate panel of Fig. 6 b shows the potential difference between CE and RE2 (E CE -E RE2 ), and the lowest panel gives the potential difference between RE1 and RE2 (E RE1 -E RE2 ). Interestingly, E CE -E RE2 (middle panel of Fig. 6 b) shows clear fluctuations during operation. Except at the beginning, these fluctuations are well aligned with the positions of the steps in E WE -E CE curve, as indicated by the derivative curve dV/dt on the top panel of Fig. 6 b. Minor fluctuations in E RE1 -E RE2 (lowest panel of Fig. 6 b) can also be observed when E WE -E CE experiences steps. These relations are labeled by green and red shaded areas for charge and discharge.
CE in this cell is a metallic Li electrode. During the charging and discharging of this cell, Li plating and stripping happen at the CE. Li plating on a metallic Li electrode happens either in the root-growing mossy Li mechanism or tip-growing dendritic Li mechanism [63] . Both the mossy and dendrite Li growth will not give rise to the regular fluctuations of voltage curve within a symmetric Li/Li cell [ 63 , 64 ]. At the beginning of Li plating, a voltage valley (middle panel of Fig. 6 b) is related to the Li nuclei formation [ 15 , 65 , 66 ]. During Li stripping off a metallic Li electrode, Li will be first electrochemically dissolved from the previously plated Li, followed by the dissolution of the bulk Li electrode. Smooth profiles can be found in the voltage within a symmetric Li/Li cell, and a single sharp peak at the middle-to-end of stripping is sometimes found due to pits formations [ 64 , 67 , 68 ]. The cell used in Fig. 6 is not a symmetric Li/Li cell, but the plating and stripping mechanism on the CE should be the same. Multiple peaks (or valleys) observed in the voltage of E CE -E RE2 and E RE1 -E RE2 in Fig. 6 b are not likely related to the growth of mossy and dendrite Li or the formation of pits as explained above. The only possible reason left is the fluctuation of Li concentration in the electrolyte 2 when adopting a porous graphitebased electrode as the WE. As explained in the first three sections, the transitions from the plateaus to slopes (steps) in the surface potential C cause wave formations and propagations in the reaction distribution C when eletrode utilization is not uniform, further provoking the waves of Li concentration in the electrolyte 2 . The steps in the output voltage at low C-rates are precisely located at the valley parts of the 2 waves due to the mass conservation of the electrolyte and charge conservation in the porous electrode. In Fig. 6 b, the positions of peaks (or valleys) E CE -E RE2 and E RE1 -E RE2 are well aligned with the steps in battery output voltage (E WE -E CE ), nicely demonstrating the behavior of 2 waves on the graphite-based porous electrode. Similar voltage fluctuations between the reference and counter electrodes in a three-electrode setup with graphite as the working electrode were previously reported [ 43 , 69 ]. The variations of Li + concentration in the electrolyte were also concluded to be responsible for these phenomena. Fig. S12 gives the simulation of E WE -E CE , E CE -E RE2 , and E RE1 -E RE2 based on Equations S16-S28 and Fig. S11. Note that the simulation of E CE -E RE2 includes some essential properties, i.e., the kinetics of the Li plating and stripping and Li + concentration in the electrolyte at the surface of CE ( 2 ), but does not include the exact activity of ions, geometry properties, microscale details of Li plating and stripping. E WE -E CE shows good agreement between experiment and simulation. E CE -E RE2 shows moderately good agreement. E RE1 -E RE2 does not show a very good agreement between the experiment and simulation because the measured difference is minor, below mV level. At the same time, the qualitative agreement is present and general fluctuations, caused by waves of 2 , are well captured by the simulations.
Once again remind, that the thermodynamic factor (OCV) and kinetic factor (reaction heterogeneity) together cause the waves of C , which further create waves of 2 and therefore peaks (valleys) in E CE -E RE2 and E RE1 -E RE2 . Generally, the reaction distribution inside the porous electrode is not uniform, which means the kinetic factor is fulfilled. The thermodynamic factor is determined mostly by the equilibrium potential curve of the electrode, and different electrode materials have different equilibrium potential curves, which influence the 2 waves. To illustrate this statement, Fig. S13 shows the potential differences among the four electrodes of EL-Cell with a porous LFP electrode as the WE under 0.1C of LFP electrode (0.16 mA/cm 2 ). The LFP electrode measured with respect to Li + /Li has a relatively flat OCV-curve except at the beginning and the end of (dis)charging. It can be seen that the E CE -E RE2 and E RE1 -E RE2 curves show smooth profiles without any fluctuations at the flat regions of the (dis)charging potential curve, indicating no waves of 2 formed at these regions. The increase of E CE -E RE2 and a valley (peak) at the end of (dis)charging are due to the changes in kinetic of Li electrode [ 64 , 69 ]. According to the model, waves of 2 can be seen near the steps in at low C-rates. This phenomenon also holds at the end of (dis)charging when voltage drops fast. Corresponding jump in E RE1 -E RE2 at the end of discharging can be seen in Fig. S12. Fig. 7 a gives the potential differences between electrodes of the graphite-based EL-Cell at various charging C-rates. Fig. 7 b gives the potential differences at two discharging C-rates. For (dis)charging at low C-rates (less than 0.2C), clear fluctuations in the potential difference of E CE -E RE2 and E RE1 -E RE2 can be observed when the battery voltage undergoes step, as shown by the green and red shade areas. For (dis)charging at middle or high C-rates, the fluctuations in E CE -E RE2 and E RE1 -E RE2 are not noticeable. The possible reason for these phenomena is the high electrolyte concentration gradient at middle or high C-rates. Such gradient causes small changes resulting from 2 waves to be invisible. In these cases, smooth profiles in E CE -E RE2 and E RE1 -E RE2 can be found when (E WE -E CE ) show flat curves for (dis)charging.
Other factors, i.e., the porosity of the electrode and separator, the separator thickness, the amount of electrolyte, etc. , will also influence the potential difference of E CE -E RE2 and E RE1 -E RE2 due to the conservation law in 2 . The roughness and neatness of metallic Li will also affect the values E CE -E RE2 . In real applications, the reaction distribution can be nonuniform in both in-depth and in-plane directions [70][71][72][73] . For the battery with large-sized electrode or high C-rates applications [ 72 , 74 ], the nonuniformity of reaction distribution will be even more pronounced at the in-plane direction, causing a more complex electrolyte concentration distribution.
Li + concentration waves in the electrolyte are investigated in the present paper using the combination of P2D-based modeling and experiments. To sum up, Fig. 8 shows a schematic illustration of the origins of the Li-ion concentration waves. As shown, the synergetic effect between the kinetics and thermodynamics causes the reaction distribution waves to mitigate the nonuniform utilization of the active material in porous electrode. The kinetic factor refers to the charge-transfer reaction heterogeneity, which causes the nonuniform utilization of porous electrodes. The thermodynamic factor refers to the OCV of the cell. The plateaus of the OCV curve tend to make the electrode utilization even more nonuniform.
In contrast, the sloping parts make the electrode utilization relatively uniform. The voltage transitions from the plateaus to the sloping parts of the OCV curve generate the waves of the reaction distribution. Such reaction distribution waves subsequently cause a local abundancy or deficiency of Li + in the electrolyte since Li-ions act as charge carriers in the electrochemical reactions. That eventually results in Li-ion concentration waves in the electrolyte. Changing thermodynamic or kinetic factors will influence the reaction distribution and eventually affect the Li + concentration waves in the electrolyte. For example, low porosities and thick electrodes will lead to larger reaction distribution waves, which cause stronger Li + concentration waves in the electrolyte.
Note that the concentration waves were not found in previous publications with P2D-based models, even though these models have been extensively used to simulate Li-ion batteries. There are two possible reasons for that. First is the model simplifications. If simplification causes inaccuracy in C then Li + concentration waves in the electrolyte can be weakened or even disappear. For example, applying uniform reaction distribution ( C ) can significantly reduce the complexity of the model [ 56 , 57 , 75 , 76 ]. However, the waves in reaction distribution and electrolyte concentration disappear, as shown in Fig. S1 in the Supporting Information. Another possible distortion is fitting the equilibrium potential or open-circuit voltage (OCV) by polynomials [ 22 , 24 , 77 ]. The resulting polynomial curves are smooth and do not contain details of the single-phase or solid-solution stages during (de)lithiation. Such a simplification eventually depresses the concentration waves in the electrolyte, as revealed in Fig. 3 c, f, and i in the manuscript.
The second reason is that concentration waves are easily overlooked when only the concentration profiles as a function of (normalized) position are shown. In most publications using P2D models, the concentration in the electrolyte is shown with respect to the (normalized) position. A limited range of discrete moments of time has been selected and compared in the same figure [ 16 , 17 , 21 , 24 , 52 , 78 , 79 ]. In this case, the selected range of moments is not enough to reveal Li concentration waves. However, if a small timestep is selected and a long simulation time is used, then some indications of Li + concentration waves can be observed. For example, there are several publications showing concentration in the electrolyte as a function of time [ 19 , 23 , 52 , 80 ]. In these cases, it can clearly be seen that the Li + concentration in the electrolyte fluctuates with respect to time at a particular position. This observation directly points to Li + concentration waves formed during operation.
The influence of Li + concentration waves on the cell performance is another important question. Two aspects can be considered here. The first one is related to a direct influence. The cell's output voltage is a summation of the equilibrium potential and overpotentials. The Li + concentration difference in the electrolyte is a critical source of the overpotential. As revealed in the present paper, Li + concentration in the electrolyte shows strong waves during operation, eventually causing fluctuations in the overpotential term [44] . The second aspect is the indirect influence, resulting from the relation between concentration waves and the reaction distribution C . C alters the utilization of the active materials inside the porous electrode, and the utilization of the active materials thermodynamically and kinetically influences C in reverse. These complex and coupled relations significantly affect the output voltage curve of the battery. For example, the phase-change indications diminish at high C-rates, as indicated by the blue shaded areas in Fig. S9a.
Finally, the findings in the present paper will also provide guidelines for porous electrode and battery design. In practical applications, achieving uniform reaction distribution is always beneficial for the energy/power output because only in this case the active materials can be fully utilized. Therefore, reaching uniform reaction distribution can be an important goal in designing high-energy/power Li-ion batteries. From thermodynamic considerations, using active material with smoothly changing OCV without flat parts resulting from phase changes helps to build a more uniform reaction distribution across porous electrodes and, consequently, maximize the energy/power output, for example, hard carbon in anodes [81] , LiNi x Co y Mn 1-x-y O 2 in cathodes [54] . Permitting some limited (5-10%) nonuniformity in utilization could be an achievable strategy for balancing the battery design and maximum extracted energy considering the kinetic factor.

Conclusions
In the present paper, Li + concentration waves in the electrolyte of a C/Li cell during the galvanostatic operation are investigated with the help of the P2D model and experimental measurements with a fourelectrode setup. The simulations reveal that Li + concentration waves in the electrolyte are generated by the waves of reaction distribution inside the porous electrode. Two factors are essential for the generation of reaction distribution waves. The reaction heterogeneity (kinetic factor) results in a non-uniform electrode utilization. The non-uniform electrode utilization tends to be mitigated when going through the voltage transitions in the OCV curve (thermodynamic factor) by forming reaction distribution waves. In the experiments, a four-electrode device (WE/RE1/RE2/CE) is used to validate the electrolyte concentration waves inside the separator region. The potential differences between CE and RE2, RE1 and RE2, show apparent fluctuations during operation, demonstrating waves in the electrolyte concentration. The alignments of the steps in battery output voltage and the fluctuations in the potential differences illustrate electrolyte concentration waves' dependency on thermodynamics and the porous electrode's reaction kinetics.
Such phenomena are not specific to graphite-based porous electrodes. It is expectable that the porous electrodes with flat and sloping parts in OCV curves (Fig. S15), i.e. spinel LiMn 2 O 4 , LiV 2 O 5 , NaCoO 2 , will also generate reaction distribution waves and electrolyte concentration waves at low-current applications. The results in the present paper reveal that one porous electrode inside a battery will influence the other (porous) electrode through the electrolyte. These results also provide an important viewpoint for understanding the reaction nature and explaining voltage artifacts in porous electrodes. In addition, the finding of the dependency of the reaction distribution on both thermodynamics and kinetics is highly relevant for the scientific community in the field of rechargeable batteries and applications. The presented results will also broaden the design prospective of porous electrodes and enhance the performance of future Li-ion batteries.

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