Modeling the kinetic behavior of the Li-RHC system for energy-hydrogen storage : ( I ) absorption

The Lithium-Boron Reactive Hydride Composite System (Li-RHC) (2 LiH + MgB$_{2}$ / 2 LiBH$_{4}$ + MgH$_{2}$) is a high-temperature hydrogen storage material suitable for energy storage applications. Herein, a comprehensive gas-solid kinetic model for hydrogenation is developed. Based on thermodynamic measurements under absorption conditions, the system's enthalpy $\Delta$H and entropy $\Delta$S are determined to amount to -34 $\pm$ 2 kJ $\cdot$ mol H$_{2}^{-1}$ and -70 $\pm$ 3 J $\cdot$ K$^{-1}$ $\cdot$ mol H$_{2}^{-1}$, respectively. Based on the thermodynamic behavior assessment, the kinetic measurements' conditions are set in the range between 325 {\deg}C and 412 {\deg}C, as well as between 15 bar and 50 bar. The kinetic analysis shows that the hydrogenation rate-limiting-step is related to a one-dimensional interface-controlled reaction with a driving-force-corrected apparent activation energy of 146 $\pm$ 3 kJ $\cdot$ mol H$_{2}^{-1}$. Applying the kinetic model, the dependence of the reaction rate constant as a function of pressure and temperature is calculated, allowing the design of optimized hydrogen/energy storage vessels via finite element method (FEM) simulations.

et al. [17]. By the combining 1 mol MgH2 with 2 mol LiBH4 a fully-reversible composite system could be obtained [22]. In this material (hereafter named Li-RHC), the theoretical absorption reaction under around 400 °C proceeds as follows [22,23] Under this assumption, the reaction proceeds in a one-step fashion, with the reactants' consumption considered to occur concomitantly. The overall absorption reaction is exothermal. It is however important to see that the decomposition of the stable MgB2 (standard enthalpy of formation of MgB2 is ∆H = -92.0 kJ•mol -1 ) is endothermal [24]. Therefore, the overall enthalpy of reaction decreases, allowing the formation of LiBH4, otherwise difficult to obtain. Thermodynamic calculations have shown that the theoretical reaction enthalpy of equation (1) amounts to -46 kJ•mol H2 -1 [17]. Absorption and desorption behaviors are, however, markedly different for this material [25], what brings even more challenges for interpreting data and greatly limits the applicability of conclusions drawn from experimental investigations and theoretical calculations to their respective case-scenarios.
This reaction is generally accepted [22,23,[25][26][27][28][29][30][31][32] as representative of the absorption process and proceeds as a one-step reaction [25,28]. However, it was first suggested by Vajo et al. [22] that between 400 and 450 °C a two-plateau region should exist. Cova et al. [23] have shown that after around 413 °C a two-plateau region can be seen. These plateaus are related to the equilibrium conditions of the LiBH4 phase and the Mg/MgH2 reaction. Although such expected transition has been observed in previous works, enthalpy and entropy values are seemingly in disagreement even for the single-plateau temperature range among different works [22,23,33]. The kinetic properties of the system are influenced fundamentally by temperature and pressure. However, they are also affected by many other factors, including, for example, the reacted fraction, use of additives (catalysts) [22,[25][26][27][28][29][30][31][32]34], microstructure, particle size distribution [35], cycle number [9], and degree of compaction [36], among others.
Envisioning the application of Li-RHC in a hydrogen storage system, a complete kinetic and thermodynamic investigation is required. With this information, it is possible to develop kinetic models that describe the reaction rate of the material as a function of the reacted fraction, the operative pressure and the temperature. One of the main challenges in the design of hydrogen storage systems based on hydrides lays on the development of numerical models that describe the phenomena occurring upon hydrogenation and dehydrogenation. Trustable models allow to evaluate the practical feasibility of the system, reducing the time and costs associated with experimental evaluations. Furthermore, the numerical development can lead to the identification of the most relevant parameters, so that the design can be optimized, either by purely numerical methods or by its combination with novel computational approaches such as machine learning [37].
This work aims to develop a comprehensive kinetic model for the absorption reaction of Li-RHC (2 MgH2 + LiBH4) with 0.05 mol TiCl3 as an additive. For this purpose, the assessment of the thermodynamic behavior is performed, as well as investigations of its kinetic properties. These results enabled the determination of the rate-limiting step of the hydrogenation process, the identification of the driving-force component and consequently, the calculation of apparent activation energy for the absorption reaction under a wide range of pressure and temperature conditions. Finally, the calculations allow the identification of a general equation that describes the kinetic behavior of the system for the chosen experimental conditions. To the best of our knowledge, this is the first time that a comprehensive kinetic model for the hydrogenation of Li-RHC is presented, which contributes with new insights for forthcoming investigations about the modeling and design of hydrogen-energy storage reservoirs. This

High-Energy Ball-Milling (HEBM)
For the preparation of the Li-RHC, LiH powder (Alfa Aesar, purity of ≥ 99.4 %) and MgB2 powder (Alfa Aesar, purity of 99 %) were mixed in a 2:1 molar ratio. Then, 0.05 mol of TiCl3 (Sigma Aldrich, purity ≥ 99.995 %) per mol of MgB2 was added. The milling process was carried out in a planetary ball-milling device (Fritsch, Pulverisette 5, Germany) using a 76 mm diameter tempered steel milling vial (66 mm high) with 10 mm tempered steel balls under argon atmosphere, a BPR (ball mass to powder mass ratio) of 10:1, for a total milling time of 20 hours (4 hours milling, followed by 1 hour wait time, repeated 4 times), at 230 RPM, with 20% of volume filling (ball volume with respect to the vial internal volume). For each milling, an amount of 5 g of powder was inserted into the milling vial. These parameters were chosen based on analyses developed and thoroughly discussed in previous work [35]. The as-milled powder present the following characteristics: particle size ranging from 10 to 70 μm, surface area of about 15 m 2 /g and a rounded-platelet-like morphology [35].
All the handling was performed under argon atmosphere and the storage of the samples was done in a continuously purified argon-filled glove box (MBraun, Germany).

Intrinsic Kinetic measurements, data handling and PCI curves
The hydrogen absorption experiments to assess intrinsic kinetic * behaviors were performed using a Sieverts-type apparatus (HERA, Canada, Canadian Patent, Serial Number 2207149 [38]) equipped with a differential pressure sensor and calibrated volumes. The internal apparatus temperature was maintained at 40 ± 1 °C at all times. The sample heating was provided by an oven surrounding the whole sample holder. The thermocouple used to assess the temperature during the experiment was located on the sample holder's outer wall.
For these measurements, around 100-150 mg of material were used to assure isothermal and isobaric conditions, as well as homogeneous concentration changes in the mass of material. With this, the mass of material can be considered as a punctual mass, avoiding the influence of heat transfer and mass transport phenomena [39]. This condition is also taken in several works [9,25,27,34,48] applying gas-solid models to fit the experimental curves, in which the used mass ranges between 100 -200 mg. The use of more mass can cause a deviation in the analysis of the intrinsic kinetic model, and mainly create a non-uniform temperature profile in the sample, so that one have neither isothermal nor homogeneous concentration changes in the used mass of material. Such condition leads to mismatch in the real intrinsic kinetic behavior, since the equilibrium pressure of the material in different parts would be different as the temperature changes, leading to a concentration profile in the material [39]. The experiments were performed starting at the 18 th cycle to ensure that the material was already stabilized in terms of capacity and kinetic behavior (see also Figure S1 and Figure  S2 of the Supplementary Material). The two sets of experiments were kinetic measurements at constant pressure (30 bar), with varying temperatures ranging from 312 °C up to 425 °C (10 measurements in regular intervals), and kinetic measurements at constant temperature (375 °C) with varying pressures, from 15 bar to 50 bar (with a step of 5 bar).
The acquired data was batch-processed with a specific Python (Python 3.7) script and further treated in OriginPro™ version 9.6.0.172 Software (Origin Lab Corporation). The fittings and statistical evaluation employed user-defined functions, which were added with the in-built tools of the program. The expression used for uncertainty propagation is shown in equation (S1) in the Supplementary Material.
The Pressure-Composition-Isotherms curves were acquired between 350 °C and 425 °C using a PCT Pro (SETARAM, Caluire, France) at the University of Pavia, Italy. A mass of approximately 150 mg was used. The amount of absorbed hydrogen was normalized as reacted fraction according to equation (2) for each of the curves individually.
where mH2, t is the mass of hydrogen absorbed at a given time t after the experiment started, mH2, max is the maximum capacity reached at the time in which the experiment was terminated.
Reacted fraction values are experimentally determined and always vary from 0 to 1.
To determine ∆H and ∆S via van't Hoff equation, the hydrogenation process of the studied material was considered as a single-step process. The start and end of the plateau region where taken, respectively, as 0.2 and 0.7. For these calculations, the mean value of the temperature Tmean throughout the experiment was considered.
The ∆H and ∆S are explicitly considered with a negative sign, since the hydrogen absorption reaction is exothermic.

Empirical Kinetic Model and General Definitions
The approach used in the present work is the Separable Variable approach as described by equation (3) [40,41]. In this model, it is possible to obtain the variation of the reacted fraction α as a function of time t by determining three variables, namely, K(T), F(P) and G(α). The name of this method of deriving a kinetic model comes from its strategy, which implies: by keeping two (undetermined) variables constant, it is possible to determine the third one. By reorganizing the variables and keeping some of those constant, all of the variables can be determined.
The G(α) term in equation (3)  (size and geometry) [42,43]. G(α) is represented by different expressions according to the gassolid reaction model and as a function of the quantity reacted fraction α [40].
These different gas-solid kinetic models belong to different categories and each implies a rate-limiting step for the overall reaction progress. In Table 1, it is possible to see a description of the gas-solid models with their names, differential form G(α), integral form g(α), acronym taken in this work and corresponding references. There is a large number of gassolid models, and thus only some of the best-fitting models considered in the present work as possible candidates are shown. For the complete table, please refer to Supplementary Material Table S1. For a more in-depth discussion of the application of these gas-solid kinetic models, please refer to the work of Puszkiel [40]. The kinetic constant can be expressed as a function of the temperature and pressure functionalities, as described in equation (4): K(T) represents the temperature-dependent term. This functionality is defined by the Arrhenius form, as shown in equation (5).
where A is the frequency factor (also called pre-exponential factor) and Ea is the apparent activation energy (of the reaction occurring in the material as a whole). The pressure-dependency component F(P) is related to the chemical reaction's driving force. This component includes the operative pressure P and the equilibrium pressure Peq. The precise form of F(P) is not easily determined and it is usually tentatively chosen among possible candidates, which are functions of parameters mentioned above. For the sake of clarity, Table 3 in Section 3.2.2 shows the tentatively used expressions in the attempt of fitting the experimental data.
To calculate the F(P) values, the Peq expression was obtained from the entropy ∆H and enthalpy ∆S values by using the van't Hoff equation as described by equation (6) in Section 3.1. In all the calculations, the temperature and pressure were taken as the mean temperature Tmean and mean pressure Pmean throughout the course of the experiment.
A simplified description of the steps to perform such calculations is provided below.
1. A set of curves under the same H2-pressure under different (but quasi-constant) T are collected. For each of these curves, K(T) and F(P) are assumed constant. By calculating the reacted fraction α from the experimentally acquired kinetic curves from 10 to 90 % of α it is possible to evaluate which reaction model G(α) best represents the kinetics of the studied system. The "reduced time method" proposed by Sharp et al. [49] and Jones et al. [50] was used for this determination. For additional information, please refer to the Supplementary Material; 2. Once a reaction model was assumed, a non-linear fit considering the integrated form g(α) is used to obtain individual sets of k(T, P) from the experimental data; 3. With these values of k(T, P), a first determination of the frequency factor A and the apparent activation energy Ea is done by considering a reaction model G(α) and F(P) as constant for each curve. The values found here are not the final ones and new values are assumed in the following steps; 4. Assuming the calculated A, Ea, and the reaction model G(α), the different pressuredependency term F(P) expressions are evaluated. After evaluation, if a suitable expression is found, new values for the frequency factor A and the apparent activation energy Ea are assumed, thus, this is an iterative process; 5. With the determination of K(T), F(P) and G(α), it is now possible to perform the data validation by plotting the calculated and the experimental reacted fraction α against time t curves.
For the determination of a general reaction model expression, some remarks should be taken into consideration. First, the models presented here were originally developed for the modeling of nucleation and growth kinetics considering liquid-solid (solidification) and solid-solid interactions [44][45][46]. Later, these same models were adapted for gas-solid interactions for hydride kinetic modeling. In the present case, all the kinetic measurement temperatures are above 270 °C, which is the melting point for the LiBH4 phase [28]. However, the fact that one of the phases is in the liquid state was not considered in this analysis. This is possible particularly because these models presuppose themselves the existence of a rate-limiting step for the reaction. The gas-solid model describes the changes of the materials upon hydrogenation and dehydrogenation based on the fitting of such models to experimental curves. An analysis of this complex system at the atomic level is out of the scope of this work. Moreover, the overall reaction, as shown in equation (1), considers that the reaction occurs as a single-step reaction. All these simplifications were taken and validated in previous works about the analysis of the rate-limiting step of this complex hydride system [9,25,48].  Figure 1 shows the hydrogenation pressure-composition isotherms measured in the range between 350 °C and 425 °C. The PCIs display a notable variation of the equilibrium pressure Peq with the amount of absorbed hydrogen (transformed fraction), which gives rise to a "sloped" plateau. The thermodynamic parameters enthalpy ∆H and entropy ∆S can be calculated from the measured equilibrium pressure Peq by applying the van't Hoff equation (6). At each temperature, the Peq was determined as a mean value between 0.2 and 0.7 of the transformed fraction α. Furthermore, a mean temperature Tmean was taken considering the whole duration of the experiment. The P0 is the thermodynamic reference pressure (considered 1 bar).

Thermodynamical Properties
By determining enthalpy ∆H and entropy ∆S values, it is possible to calculate an equilibrium pressure Peq as a function of the temperature T.
A comparative evaluation of the results for ∆H and ∆S and the pressure ranges reported in the literature is seen in Figure 2. The data points in Figure 2 were taken from the works of Vajo et al. [22], Cova et al. [23], Puszkiel [33]  instead being shown. Still, the differences between the recalculated values and the reported values were significantly different only for the single plateau for the work of Cova et al. [23]. A complete table with the used pressure values, the fitting parameters, the calculated and reported values is available in Table S2 in the Supplementary Material.
The straight lines represent the linear regression performed with the points shown in the diagram. The dashed line at around 412 °C is an approximation for the temperature range after which the reaction would not occur as a one-step reaction. Vajo et al. [22] has reported that the ∆H and ∆S for this system amounts to -40.5 kJ•mol H2 -1 and -81.3 J•K -1 •mol H2 -1 . Moreover, it has been proposed that the rise of a two-plateau reaction could possibly be observed for higher temperatures [22]. The data provided for the 450 °C PCI on that work could not resolve this issue. Later on, works from Puszkiel [33] and Cova et al. [23] have both been able to measure the presence of two plateaus for temperatures over 400 °C to 412 °C, respectively. In this range of temperatures, the system presents a region between the lower (LiBH4) and upper plateau (Mg/MgH2) in which it is theoretically possible to have coexistence of LiBH4 and Mg in equilibrium conditions. This twoplateau behavior, however, is not clearly visible in the results presented in this work. This fact can be attributed to the differences in the equilibrium conditions resulting from differences in starting materials, handling, processing, among other experimental conditions. It should be expected that such a two-region plateau exists even in our material, but it may be possible that this transition occurs at higher temperatures.
As seen in Figure 3.a), for the hydrogenation process, the obtained values for ∆H and ∆S in this work are -34 ± 2 kJ•mol H2 -1 and -70 ± 3 J•K -1 •mol H2 -1 respectively. The obtained fitting goodness (R²) of 0.994 shows a proper correlation. It is important to notice that the redmarked point at 425 ºC is not considered for the linear fitting. These values are in good agreement with the results published by Vajo et al. [22] for the temperature ranges below 400 °C. In Figure 1, for the curve at 425 °C, the "bump" seen around 0.55 and the steep increase in the pressure with increased reacted fraction suggests the presence of a second plateau. However, the difference observed here is not comparable to the change that has been reported in the work of Puszkiel [33] and Cova et al. [23]. Additionally, it should be noted that the point at 425 ºC (Figure 3.a)) presents a positive deviation over the fitted curve, in agreement with the work of Vajo et al., where the positive deviation is seen at 450 °C [22]. This positive deviation can be ascribed to the presence of a two-plateau region. This issue is not still clear and would require further investigations that are beyond the scope of this work.
In the present work, for the sake of preciseness, a conservative assumption is made to guarantee the validity of the kinetic model in a range in which the mutual hydrogenation of MgB2 and LiH to LiBH4 and MgH2 occurs. The temperature range of the kinetic data is limited to 412 °C. Above 412 ºC it may be possible that the system undergoes a different reaction pathway, and the kinetic analysis would not be representative. Therefore, the 425 °C kinetic results are, from now on, excluded from the kinetic modeling dataset.
With the calculated values of ∆H and ∆S it is possible to draw a pressure vs. temperature diagram (an exponential form of van't Hoff equation) as seen in Figure 3.b). In this diagram, it is possible to see the equilibrium conditions for the Li-RHC system considering only the results found in the present study. Additionally, the points considered for the development of the kinetic model are shown as round markings. Additionally, the X-markings seen in Figure 3.b) are the experimentally measured mean plateau pressure values. It is possible to see that there are only small deviations between the values obtained from the van't Hoff equation and the experimental results used for the fit.
The general expression for calculating the equilibrium pressure Peq is given in where R is the ideal gas constant, and Tmean is the mean temperature throughout the experiment.

Determination of G(α): Gas-Solid model
Different reaction models have been developed to represent hydrogenation and dehydrogenation reactions. The models here considered for this purpose can be seen in Table  S1 of the Supplementary Material. To deduce which of the reaction models best describes our system under the studied experimental conditions, the general approach is to start evaluating which of the integral equations g(α) (see Table 1) best correlates with the obtained experimental data. In order to do so, the so-called "Sharp and Jones" method (named after the original works of Sharp et al. [49] and Jones et al. [50] and also known as "Reduced Time Method") has been shown to be the most efficient tool, as it gives more parameters for evaluation that help determining how good is the agreement of the experimental data to each of the reaction models. A general description of the employed method has been given in the work of Puszkiel [40].
The evaluation of the best-fitting model is based on the values of three parameters; the coefficient of determination R² closest to 1, a Y-axis intercept closest to 0 and a slope closest to 1 [40,49,50]. All results are summarized in Table 2. Further information on how the method was implemented can be seen in the Supplementary Material, along with the figures for all the fittings done in this work ( Figure S3).
All the fittings were performed for the measured temperatures within the range of 312 °C to 412 °C under 30 bar of initial H2-pressure for each of the models. A graphical summary of the results in regards to the models in each temperature can be seen in Figure S4 of the Supplementary Material.
The obtained results indicate that there are 5 different models that rank sufficiently well in regard to R² values, namely, JMAEK with n = 1, JMAEK with n = 1.5, CV, 2-D diffusion limited and 3-D diffusion limited models. However, as slope and Y-intercept agreement are taken into account, both diffusion-limited models cannot be considered as suitable candidates. Regarding the three remaining models, JMAEK with n = 1.5 has insufficient proximity of the target values for slope and Y-intercept. Slightly closer to the targets is the CV model. Still, the JMAEK with n = 1 has, comparatively, the best match for the two parameters, especially considering Y-intercept.  Taking into account that the JMAEK with n = 1 model presents the highest R² values for all the temperatures (range from 0.996 to 0.999) and ranks rather well in the two other parameters, it is from now on assumed to be the most suitable model to describe the kinetics of the system.
For pristine Li-RHC, different authors reported differing models and even rate-limiting steps. In one of the first works related to the determination of rate-limiting steps and reaction models, Wan et al. claimed that the rate-limiting step should be the diffusion of species through the product layer [51]. Differently, Bösenberg et al. argued that an interface-limited model would be a better representation of the system [27]. Later on, Puszkiel et al., in two different studies, argued in favor of an interface-controlled rate-limiting step [34,48]. Le et al. also reported an interface-controlled reaction [25]. Two models with interface-controlled kinetics were presented in these studies; the three-dimensional contracting-volume interfacecontrolled (3D CV) [27,34] and one-dimensional growth with interface-controlled reaction (JMAEK with n = 1) [25,48]. Although different models have been proposed, it should be noted that these models imply the same rate-limiting step. Furthermore, comparing these previous studies, some remarks should be taken into consideration. First, in some earlier works the JMAEK with n = 1 was not taken into account. Second, the method for determining and evaluating the fitting goodness of the kinetic models became more robust with time, since it became increasingly common practice to employ the Sharp and Jones method [40,[48][49][50].
Concerning the Li-RHC with additives, different models were proposed for the system. Puszkiel et al. reported that, with the addition of 1 % mol of TiO2, the model that best represents the system is the interface-controlled one-dimensional growth (JMAEK with n = 1) [48]. Le et al. reported that for the Li-RHC with the addition of 0.00625 mol of (3TiCl3·AlCl3), the bestdescribing model changes in relation to the pristine material and becomes the 3D contractingvolume interface-controlled [25].
Particularly diffusion-controlled) and the interface-controlled three-dimensional contracting-volume (3D CV interface-controlled) [27]. In their work, though, the JMAEK with n = 1 was not considered as a candidate. Additionally, the so-called reduced-time method based on the works of Sharp et al. [49] and Jones et al. [50] was not employed at the time.
In more recent studies for the 0.05 mol TiCl3-added Li-RHC system it was found that the JMAEK with n = 1 best represents the experimental data [9]. This determination of g(α) was made based on different temperatures between 325 °C and 425 °C under 30 bar of initial H2-pressure [9].
In the present work, by measuring the absorption kinetics in smaller steps (total of 10 measured curves between 312 °C and 425 °C) and by using the results after the 18 th absorption cycle, it is reasonable to assume that both the influences of the calculation errors and of the change in hydrogen capacity with the first (de)hydrogenation cycles (altering the time to reach 90% of the reacted fraction) have been accounted for or reduced, respectively.
Considering that the JMAEK with n = 1 model describes a one-dimensional interfacecontrolled reaction, the results reported on [25,27,34,48] are in line with what is being proposed in this work.
The experimental results for reacted fraction α against time are shown in Figure 4.a) and Figure 4.b). The solid lines represent experimental kinetic data at different temperatures under 30 bar of initial H2-pressure and for different pressures at 375 °C, respectively. The nonlinear fit was performed to determine the k-values for the integrated JMAEK expression with n = 1 [44][45][46], as shown in equation (8): The non-linear fitted curves are presented in Figure 4.a) and b) as dashed lines. Both graphs indicate clearly that higher pressures and higher temperaturesas usually expectedresult in faster kinetics. the fitting does not produce curves with inflections, which is observed clearly in the experimental curves for 325 °C and 337 °C in Figure 4.a). Such an aspect has been previously identified also in other works performed on this material [9]. The fitting parameters can be seen in Table S3 and S4 of the Supplementary Material. A different kind of misfit between experimental data and fitted results is seen for the kinetic curves at 375 °C. For 15 and 20 bar pressures, the non-linear fits seemingly underestimate the hydrogen uptake during the first minutes of reaction. However, the increase in the reacted fraction for the experimental curves is likely related to an experimental artifact due to the time necessary to close the valve that initially sets the pressure difference between reference and sample holder volumes to zero. For a more in-depth description of the internal workings of the measurement apparatus, see [38]. It should be considered that the outline of the fitted curve presents a more realistic representation of the kinetics of the material under these experimental conditions than the experimental curve itself. Considering all the described effects, the deviations between the model and the experimental results are negligible. Therefore, it is possible to conclude that the proposed gas-solid model (equation 8) is in good agreement with the experimental results.

Determination of temperature K(T) and pressure F(P) functionalities
With the values from the k(T, P)-semi-empirical kinetic constants obtained from the non-linear fits of the experimental results from different temperatures under 30 bar initial H2-pressure, the ln k against the inverse mean temperature are plotted in Figure 5. here found for the frequency factor A and the apparent activation energy Ea were respectively 100 ± 2 kJ•mol H2 -1 and (8.1 ± 2.6)·10 4 s -1 .
In order to take into consideration the influence of the experiment pressure P and the mean equilibrium pressure Peq at each temperature, the two found parameters, i.e., Ea and A, are still to be corrected by the driving force component, F(P).
For the evaluation of the different models for F(P), equation (5) is divided into both sides by F(P) and combined with equation (4) resulting in Now, by applying natural logarithms to both sides in (9) With this linearized equation, it is possible to determine a frequency factor A and an apparent activation energy Ea that take the pressure dependency into account by utilizing the values for k(T, P) (previously obtained from the fits shown in Figure 4), and each of the F(P) expressions. For the F(P) term, different functionalities are assumed, deducted either from investigations about the thermodynamic behavior or different tentatively proposed ratios and relations between P and Peq. The expressions used in the present work are presented in Table  3.

Limited
Category/Conditions Expression Name Reference To take the F(P) component into account, each expression is considered individually and a new plot, analogous to Figure 5, is built, along with each linear fit. This yields new values for Ea and A for each of the F(P) functionalities, along with the R² values for each fitted curve.
Each new set of values for Ea and A are individually considered (see also Supplementary Material Figure S5). By using the k-values obtained from the fittings shown in Figure 4 (see also Figure S6 and Figure can be obtained, for which the aforementioned values are applied. Here, each F(P) expression (along with its results A and Ea) can be checked for its validity by plotting and interpreting it graphically. By analogy with a linear equation, the equation's left-hand side is taken as a dependent variable and F(P) as an independent variable. It follows that the best correlation would necessarily have a Y-intercept as near as possible to 0, a slope nearest to 1 and a determination coefficient R² closest to 1. For this particular fitting procedure, the linear coefficient is set to 0.
The results of this procedure for each F(P) expression can be seen in the Supplementary Material in Figure S6. The summarized results of these fitting parameters can be seen also in the Supplementary Material, in Figure S7. Comparatively, two F(P) expressions rank for R² much higher than the others, namely, F2 and F3, with R² values of around 0.990. These two best-ranking F(P) expressions are shown in Figure 6.a) and b). The slope for F2 and F3 is respectively 1.05 ± 0.03 and 1.04 ± 0.03. These values are considered sufficiently close to 1, not only because of the method's uncertainty, but also for the comparative evaluation of the results.
However, just by comparing the parameters of the fitting, it is not possible to distinguish which one of the two expressions is the best choice for F(P).
By visual evaluation, it is possible to see that for F3 the deviations of the linear fit are more likely to be spread symmetrically around the fit curve. For F2, most of the deviations are for the high-pressure points and they all have a positive deviation.
As this deviation is more systematic in F2 than in F3, F3 is being favored, as it is less likely to present substantial deviations when calculated for their α values. Another aspect to take into account are the physical phenomena that were considered to propose said driving force expressions.
However, only limited information is currently available in the literature of why these expressions can fit the data of different pressures, as shown in our results. The driving force term of a chemical reaction is usually based on the thermodynamic activity. Considering the hydrogen gas as ideal, for a gas-solid reaction in a metal hydride, the thermodynamic activity is defined as the applied hydrogen partial pressure P divided by the standard pressure P0 [69]. However, already for simple gas-solid hydride forming reactions, a wide variety of pressure dependence relations are used [58,59,68]. This fact is related to the different hydride forming materials, experimental conditions and rate-limiting steps associated to the F(P) [40].In the case of complex hydrides, it is well accepted that in general the relation between the applied pressure P and the equilibrium pressure Peq can describe the driving force term for the formation of a metal hydride. For instance, in the case of NaAlH4, the functionality of the pressure was empirically determined, and the best function was determined as the first order Taylor series' approximation (centered around Peq) of the change of the free energy for the hydrogenation process, i.e., (P -Peq)/Peq [57]. Therefore, in the herein investigated and rather complex 2 LiBH4 + MgH2 hydride system, the nature of the hydride forming materials, experimental conditions (set up) and the determined interface-rate limiting step leads to the F(P)= (P -Peq)/Peq as a best-fitting expression based on the above mentioned Taylor series' approximation. Due to its complexity, futher analysis regarding the physical meaning of the driving force term is beyond this work's scope.

General Expression and Data Validation
Considering the calculated apparent activation energy Ea, the frequency factor A, the driving force expression F3, i.e., F(P) = (P -Peq)/Peq, and the reaction model as JMAEK with n = 1, it is possible to propose a general expression for the studied system within the stated limits of temperature and pressure, i.e. from 325 °C to 412 ºC and from 15 bar to 50 bar. Starting from equation (3), substituting its components with K(T) as presented in equation (5), introducing the suitable expression for JMAEK with n = 1 (as seen in Table 1) and assuming F3 as the driving force expression, equation (12) is obtained, which is the differential form of the general kinetics expression.
In this expression, the reacted fraction α varies between 0 and 1, the frequency factor A is given in (s -1 ), the apparent activation energy Ea is given in (J•mol H2 -1 ), the temperature T is given in (K), the pressure P is given in (bar) and the time t is given in (s).
By plotting equation (13) with the transformed fraction α against the time t, it is now possible to compare the calculated and the experimental kinetic curves. Figure 7.a) shows the validation plot for the experiments with different temperatures (in the 325 to 412 °C range) under 30 bar of initial H2-pressure. In it, the solid lines represent the experimental data and the dashed lines, the calculated values. The insert on the bottom-right corner displays the same data but is limited to the higher temperature curves and the first 60 minutes of reaction.
For most of the temperature range, a very high degree of correlation between calculated and experimental results is achieved throughout the reaction. However, for the lowest temperatures investigated, e.g., 325 °C and 337 °C, the experimental curves present a perceptible inflection in the first hours, which is related to the significant complexity of the system reaction mechanism. This reaction modeling cannot capture this behavior, as the equations that are used imply a monotonic behavior and a curve that is always concave. Concerning these curves, it is visible that the kinetic reaction rate (curve slope) of the model is higher on the first minutes of reaction, slightly overestimating it when compared to the experimental data. After the inflection point, for most of the temperatures analyzed, the calculated reaction rates now are slightly underestimated. Before reaching the saturation point, the two curves again agree to a significant extent.
The coefficient of determination (R²) for the calculated model and the experimental results was calculated according to equation (S2) of the Supplementary Material (from 0 to 0.99 reacted fraction). The values give a more objective frame of comparison between the curves' goodness-of-fit and possibly serve as a benchmark for other works for the future irrespective of the material used or reaction model proposed. Particularly for the absorption reactions at 375 °C all the R² values were above 0.996.
Although small visible differences between the model (equation 13) and the experimental results can be seen, it should still be considered that the applied model has been successful in describing the kinetics of the system, as these deviations can be considered minor in the frame of the course of the reaction. In Figure 7, it is possible to observe that the model shows a quite good agreement between the calculated and the experimental results with fitting goodness ranging between 0.97 and 0.99.

Thermodynamic stability and kinetic behavior: Isokinetic contour graphs
An additional interpretation ofthe thermodynamic and kinetic results presented here is to draw isokinetic contour lines in a pressure-temperature diagram, as shown in Figure 8. The drawing of the isokinetic contour lines helps to understand how system conditions change influences the kinetic constant. For this purpose, a convenient approach is to represent the expression of the contour lines for different values of P as a function of T with chosen values of the kinetic constant k.
By rearranging equation (4), combining it with equation (5) as a function of temperature, making pressure the dependent variable, and implementing the F3 expression, it is possible to obtain the expression shown in equation (14). For plotting the isokinetic lines, k-values are arbitrarily chosen. Here, A is the frequency factor (1.8 · 10 8 s -1 ), the apparent activation energy Ea (146 kJ•mol H2 -1 ) and R is the ideal gas constant. The equilibrium pressure Peq is calculated by equation (7).
This expression yields the contour lines seen in the inset graph of Figure 8 and in Figure S8 of the Supplementary Material. In both, the hatched region (below equilibrium curve) represents the equilibrium conditions that favor the stability of LiH + MgB2, and the region above, in which the stable phases are LiBH4 and MgH2. This equilibrium line has been calculated using the van't Hoff equation (equation (7)) with the values of ΔH and ΔS obtained in the present work.
The curves drawn in the absorbed state region (LiBH4 and MgH2) are the isokinetic lines obtained from Equation (15). The values in Figure 8 were chosen to identify the experimental conditions in which at least two points are nearly intercepted by a single isokinetic line. In Figure S8, the values were chosen in order to show how the kinetic constant k(T, P) varies as a function of T and P. The region in which the model is expected to reliably describe the kinetic behavior of the Li-RHC under absorption conditions is schematically shown in Figure 8 as a square limited by the temperature and pressure conditions for which this model is deemed valid. Thus, inside this region, the calculated isokinetic contour lines are solid, and outside the region, they are dashed lines, indicating that, in principle, the shape of these curves can be known with reasonable precision only inside this validity region.
The lower temperature border is around 325 °C, since, as shown in Table 2 (and in Figure S3 and Figure S4 of the Supplementary Material), the JMAEK with n = 1 does not fit well the kinetic curve at 312 °C. On the higher-temperature side (above 412 °C), the temperature range is limited by the change of the reaction pathway, resulting in different equilibrium conditions. In relation to the pressure limits, it is possible to assume that in the lower pressure range, if enough separation between the system conditions and the equilibrium condition (Peq curve) occurs, the model can still represent the kinetics of the system consistently. However, no experimental validation has been performed below 15 bar, not only because the times for the kinetics would be significantly large, but also because one of the goals of this work is to describe with enough detail the region (for conditions of P and T) that is of interest for potential applications of this material in energy storage systems. In the highpressure region (above 50 bar), it is expected that the results obtained from the model can still represent the experimental results to some extent. However, as kinetics get faster in experimental conditions with a more significant driving force, it becomes increasingly harder to evaluate experimental results, as relatively small experimental uncertainties and experimental artifacts of measurement (valves opening/closing times and pressure transducer stabilization, inter alia) can lead to a significant change in the outcome and alter the interpretation of the results. For each of the isokinetic curves drawn in the inset graph in Figure 8, the nearly intercepted experimental curves under these conditions are shown in the main graph. As these matching curves show very similar curve outlines, it has been shown that the found isokinetic line expression is able to reproduce the kinetic behavior of the material inside the aforementioned validity region quite well.
Still, it should be taken into consideration that the pressure and temperature ranges that delimit the validity of the model in this work are conservatively taken. That is, while the physical explanations provided to the kinetic model should hold true only in these ranges, extrapolation of the model to close-enough neighboring regions still yield excellent prediction capabilities (not shown here), showing the mathematical robustness of the model even for extended temperature and pressure conditions.

Conclusion and Outlook
In this work, both thermodynamic and kinetic experimental data are analyzed in detail for the LiBH4/MgH2 Reactive Hydride Composite (Li-RHC) with 0.05 TiCl3 additive under absorption conditions. The obtained results in the range between 350 °C and 400 °C for ∆H and ∆S are -34 ± 2 kJ•mol H2 -1 and of -70 ± 3 J•K -1 •mol H2 -1 , respectively, in good agreement with previously reported values.
To the best of our knowledge, it is the first time that a comprehensive kinetic model under absorption conditions is developed for the Li-RHC by using the separable variable method. Applying this model, the effects of temperature, pressure and transformation of the hydride forming material are considered in the rate expression. The results indicate that the transformation of the forming hydride is limited by the movement of the not-hydrogenated/hydrogenated material interface, described by the one-dimensional interface-controlled model with a fixed number of nuclei and constant interface velocity, also known as JMAEK with n = 1 model. After taking the driving force component into account, the apparent activation energy Ea and pre-exponential factor A are, respectively, (1.8 ± 1.0) · 10 8 s -1 and 146 ± 3 kJ•mol H2 -1 .
The developed model successfully describes the intrinsic kinetic behavior of the system, with only minor observed deviations. For the considered range of temperatures and pressures, the model shows fitting goodness ranging from 0.97 to 0.99. On-going research is being done to comprehensively describe the desorption kinetic behavior of the here studied material. Moreover, the developed models will be applied to model pilot-plant sized energyhydrogen storage tanks through FEM simulations to better understand and optimize their designs.
The model presented here indicates that absorption times for the studied material are in the range of hours. While faster kinetics is a desirable trait for hydrogen storage materials, its influence on the performance of hydrogen storage systems depends heavily on the size and on the design choices for such systems [39]. For instance, intrinsic kinetics becomes increasingly unimportant as the size of the system increases. Still, the kinetic behavior is a core element of the functionally of these systems and needs to be properly described in order to obtain accurate descriptions of the combined effects that occur in these bigger-scale systems. However, faster kinetics have been demonstrated for Li-RHC in previous works. It has been previously demonstrated, that by producing ball-milled powders with some additives, it was possible to obtain nanostructures with in-situ-formed catalysts such as LixTiO2 and AlTi compounds that enable the completion of the hydrogenation reaction with times of less than 30 minutes under 400 °C and 50 bar [25,48]. Our expectations are, that since in these cases the rate-limiting step has been also identified as interface-controlled, the modelling approach here presented can be similarly applied to describe these materials with improved kinetic behavior.

Acknowledgements
Table S1 -Gas-solid reaction models evaluated, with its description, differential form, integrated form, acronym used in this work and references for its description 1 and/or use.  The propagation of error was handled by applying 6 7 in which is a function of variables (e.g., 1 , 2 , 3 …), is the best error estimate of and 9 is the best error estimate of the calculated value. It was assumed that the uncertainties of 10 the variables are all independent of each other. Unless stated otherwise, the statistical 11 uncertainties and error bars shown represent the standard error of the measured or calculated 12 value.

14
Coefficient of Determination (R²) Calculation 15 16 The coefficient of determination (R²) has been used in this work to assess how well the fitted 17 or calculated curves could adjust the experimental results. For the fittings performed, the 18 standard formulation of the software Origin Pro™ was used. For the calculation between the 19 experimental data and the curves calculated from the model, equation (S2) was used 20 21 in which i is the index, N is the number of measurement points, y exp is the experimental value 23 for reacted fraction and y calc is the calculated value for the reacted fraction. 24 To determine the R² of each of the curves in the different experimental conditions in which they 25 were obtained, the reacted fraction between 0 and 0.99 was used.   Vajo et al. [22] Single 603,15 0,00166 5,50 Vajo et al. [22] Single 636,15 0,00157 8,50 Vajo et al. [22] Single 673,15 0,00149 12,90 Vajo et al.

Additional Information on Methods: Sharp and Jones method implementation 47
The so-called "Sharp and Jones method" is a mathematical data treatment strategy in 48 order to ease the identification of the isothermal reaction kinetic model which would best 49 describe experimental kinetic data. Such method has been popularized by the works of Sharp 50 et. al. [48] and Jones et. al. [49] and is also known as "Reduced Time Method". Initially used 51 for solid-state reaction kinetic model identification, it can also be applied also for gas-solid 52 reaction kinetic models. The main advantage of such strategy is that it gives two additional 53 parameters which can eventually help the evaluation of the results and improve the chances 54 of correct selection of a suitable kinetic model to represent experimental data. 55 First, consider the integrated g(α) = k • t expression ( Table 1 of the main article file), in 56 which α is the reacted fraction, k is a kinetic constant for a given temperature-pressure pair, t 57 is the time and t0.5 is the time a reaction takes to reach 0.5 of reacted fraction α. As the kinetic 58 constant k does not change with the variation of t or α, it is possible to obtain equation (S3), 59 which is independent of the value of k. As each reaction model has a determined g(α) 60 expression, one can easily find the value of g(0.5) and [t0.5]theoretical.