Parameter Identiﬁcation of a Quasi-3D PEM Fuel Cell Model by Numerical Optimization

: Polymer electrolyte membrane fuel cells (PEMFCs) supplied with green hydrogen from renewable sources are a promising technology for carbon dioxide-free energy conversion. Many mathematical models to describe and understand the internal processes have been developed to design more powerful and efﬁcient PEMFCs. Parameterizing such models is challenging, but indis-pensable to predict the species transport and electrochemical conversion accurately. Many material parameters are unknown, or the measurement methods required to determine their values are expensive, time-consuming, and destructive. This work shows the parameterization of a quasi-3D PEMFC model using measurements from a stack test stand and numerical optimization algorithms. Differential evolution and the Nelder–Mead simplex algorithm were used to optimize eight material parameters of the membrane, cathode catalyst layer (CCL), and gas diffusion layer (GDL). Measurements with different operating temperatures and gas inlet pressures were available for optimization and validation. Due to the low operating temperature of the stack, special attention was paid to the temperature dependent terms in the governing equations. Simulations with optimized parameters predicted the steady-state and transient behavior of the stack well. Therefore, valuable data for the characterization of the membrane, the CCL and GDL was created that can be used for more detailed CFD simulations in the future.


Introduction
Ongoing climate change and the need to reduce greenhouse gas emissions require new technologies to meet global energy needs with renewable sources. Polymer electrolyte membrane fuel cells (PEMFCs) powered by green hydrogen show a performance that is well suited for automotive applications due to their high power density, high efficiency, and high startup availability [1]. However, due to the low operating temperatures, water management is a challenge. On the one hand, sufficient hydration of the membrane must be ensured for proper conduction of hydrogen ions, and on the other hand, flooding of the microporous layers and GDLs must be avoided.
Mathematical models have been developed to better understand the internal processes, such as the formation of liquid water and its transport, and to design more efficient and powerful fuel cells (FCs). In addition, predictive and reliable models are used to forecast the internal conditions and optimize their operation strategy, e.g., optimization of purging strategy due to nitrogen and water accumulation in the anode [2]. Therefore, the parameterization of the models is crucial. The processes and chemical reactions inside a fuel cell depend on many physical phenomena such as mass/heat/charge transport and electrochemistry. The available models can be categorized in [3,4]: (1) the degree of physical modeling (analytical, (semi-)empirical, mechanistic) and (2) the spatial resolution of the models (0D, 1D, 2D, quasi 3D, 3D). Springer [5] developed one of the first 1D semi-empirical models for that the multi-verse optimizer provided better and faster results than other applied optimization algorithms, such as GA, bee colony algorithm, or hybrid adaptive differential evolution (HADE). Mo et al. [36] applied a hybrid differential evolution algorithm to optimize the parameters of an electrochemical PEMFC model. Measurements of the fuel cell stack U-I curve and species inlet conditions were used as input parameters for the optimization. Similar, Sun et al. [37] and Gong et al. [38] used the HADE algorithm with a novel scaling factor and improved dynamic crossover probability for the parameterization of an electrochemical PEMFC model. The bee colony foraging mechanism enhances the local search. The results showed better identification of nonlinear parameters than other applied methods. In [39] they proposed a dynamic differential evolution algorithm with collective guidance factor to find the parameters of a mechanistic model.
This work aims to demonstrate the parameterization of a hybrid 3D analytic-numerical PEM fuel cell model with a standard measurement procedure for a fuel cell stack. A differential-evolution optimization algorithm and Nelder-Mead simplex algorithm from the open-source Python package scipy.optimize [40] were applied to find the best set of parameters in the bounds of known physical values. Six measurement series-three load steps and three polarization curves-with different operating temperatures and gas inlet pressures were available for optimization and validation. It was the specific challenge of the discussed model setup that the fuel cell stack was operated at low temperatures of 52-55°C (cooling inlet). The fuel cell stack system is self-humidifying as there is no external humidifier. The low operating temperatures were chosen to ensure a stable operation of the stack, because temperatures above 60°C would unnecessarily dry out the membrane and could damage it. The measurements showed strong temperature dependence, even with only 3 K temperature difference at the cooling inlet, a significant difference in the polarization curves was visible. Therefore, special attention was paid to the temperaturedependent terms in the governing equations. The values obtained for the parameters will be used to parameterize a more detailed 3D CFD model in the future.

Numerical Model and Governing Equations
The transient, single-phase, isothermal, and hybrid analytical-numerical 3D model of a PEMFC implemented in AVL CRUISE™ M software was used for the simulations. The model, developed by Tavčar and Katrašnik [13], assumes the gas species in the channel flow direction (z-direction) to be 1D and solved numerically and the species concentrations perpendicular to the gas flow in the channel are calculated analytically in 2D (x, y-direction). The pressure drop in the GDL perpendicular to the gas flow in the channels is neglected. The model for the gas flow is coupled with an electrochemical sub-model. The electrode kinetics are reduced to the oxygen reduction reaction (ORR) of the cathode catalyst layer (CCL), which is assumed to be infinitesimally thin. The schematic overview of the modeled fuel cell is shown in Figure 1. A representative unit of assumed parallel channels and ribs is extracted. In the direction of the gas flow in the channel, this unit is discretized into an arbitrary number of slices with length L. The domains of the slice are as follows: cathode and anode channel, GDL parts under-channel, GDL parts under-rib, thin CCL and membrane. While the original model of Tavčar and Katrašnik [13] neglected the electrical conductivity of the GDL, the current implementation considers its ohmic losses. Liquid water is assumed to be present only inside the membrane. Its transport by water diffusion and electro-osmotic drag is assumed as being perpendicular to the membrane plane. Convection due to pressure differences between anode and cathode and a gas crossover through the membrane were not considered. The water content of the membrane is in equilibrium with the relative humidity of the gases mixture at the interface with the CL/GDL. A more detailed view on the assumptions made during the development of the numerical model can be found in [13].

Governing Equations
The full set of equations of the hybrid 3D numerical model can be found in [13]. The following equations are highlighted because some of their parameters are used as design variables for the numerical optimization.

Species Transport in Gas Channels
The gas is assumed to be ideal and the species transport is modeled by Stefan-Maxwell diffusion. Only water vapor is modeled; liquid water is not considered in channels and GDLs, resulting in a single-phase approach. The species transport equation is divided into a 2D and a 1D part. In a channel slice it reads as follows [13]: With the concentration of species c in the channels and the velocity v in channel direction (z-direction). V represents the potential flow part of the gas flow in xand ydirection. Due to steady-state assumption of the gas dynamics, the time derivative of the species concentration in Equation (1) is set to 0. The binary diffusion coefficient D and the average species concentrationc are constant in one slice. As there is no liquid phase, the relative humidity can reach values higher than 1. This was an indirect indication of liquid water formation.

Species Transport in GDLs
No gas flow (zero velocity) is assumed in z-direction in GDLs. Therefore, Equation (1) simplifies to dc(x, y) where D GDL is the diffusion coefficient in the GDL. It is derived from the binary diffusion coefficient in gas channels by dividing D through the tortuosity τ GDL of the GDL

Electrochemical Reactions
The Butler-Volmer equation (BVE) is a relationship between the potential of an electrochemical system and its current density. The BVE states that the current density increases exponentially with activation overvoltage. In an electrochemical system such as a fuel cell, two independent BVEs for ACL and CCL can be employed. As the electrokinetics are reduced to the CCL by the model used, only the following BVE was applied [13,41]: where i 0 represents the exchange current density, p O 2 and p H 2 O are partial pressures of oxygen and water, λ stoichiometry ratio, α transfer coefficient, n number of exchanged electrons, F Faraday's constant, T temperature and η act activation overpotential of the electrode. The exchange current density i 0 is related to the reference exchange current density i 0,re f at reference temperature T re f and pressure p re f [13,41] With p O 2 and p H 2 O partial pressures of oxygen and water, b O 2 and b H 2 O oxygen and water exponent and E act,i 0 the activation energy of exchange current density.

Ionic Conductivity
The ionic conductivity of the membrane was calculated using the following Arrhenius approach with a reference ionic conductivity σ ion,re f at reference temperature T re f and reference membrane water content C w,re f [5]: With C w actual water content of the membrane, E act,σ ion activation energy, and T mem actual temperature of the membrane.

Numerical Optimization
Optimization algorithms search in general for the minimum of an objective function f (x) with one or more design variables x from the feasible design space X . [42] With x = [x 1 , x 2 , . . . , x n ] for n-dimensional design variables. x can be varied to find the minimum of the objective function f (x).
There are two general optimization strategies, local optimization and global optimization. The local optimization aims to find the minimum near a starting point and the solution must not necessarily be a global minimum. In general, mathematical models of fuel cells are nonlinear, numerically stiff and non-differentiable. Therefore, optimization algorithms that rely on the derivative of the objective function are not applicable. Direct methods (e.g., Powell's Method and Nelder-Mead Simplex Algorithm (NMSA)), stochastic methods (e.g., simulated annealing, natural evolution strategies) or population methods (e.g., genetic algorithms, differential-evolution, particle swarm optimization) are applied in case of non-differentiable objective functions. Whereas direct methods do not necessarily find the global minimum, stochastic and population methods search the given design space for the global optimum. Stochastic methods rely on pseudorandom strategies, such as Latin Hypercube Sampling (LHS) or Sobol sampling, to sample the design space and to escape local optima to enhance the probability of finding the global optimum. In contrast to direct and stochastic methods, population strategies do not move a single point toward the optimum, but use many design points (individuals) distributed throughout the design space. The information of the distributed design points is exchanged between the points to find the global minimum [42,43].
The scipy.optimization module from Python provides several local and global optimization algorithms. The differential evolution algorithm [43] gave the best and fastest results as a first indication of where to find the global minimum. The best combination of design variables from the differential-evolution optimizer was used as a starting point for a local optimizer, the Nelder-Mead Simplex Algorithm [44], to find the optimum more accurately.
As the simulations are transient, the objective value Y Obj is the integrated sum of the difference between the measured and the simulated voltage or current. As shown in Section 2.3, each of the three available measured polarization curves is based on six measured operating points. In the simulations, a current ramp from 0.02-1 A/cm 2 with a total time of 720 s (which corresponds to a slope of 1.4 mA/cm 2 per second) was defined. The difference of the simulated and measured voltage was then integrated over time to get the objective variable. For U meas (t) between the measured operating points, piecewise linear interpolation was applied.
The second set of available measured data were three load steps with voltage driven stacks, thus the difference between simulated and measured current was integrated.
Y Obj,LS = |(I meas (t) − I sim (t))|dt (9) All available measurements were used to optimize the parameters. Therefore, an appropriate combination of the objective value with weighting factors needed to be implemented. The results from the load steps were weighted with 2/3, whereas the polarization curves were weighted with 1/3. The three load steps and three polarization curves were weighted equally among themselves. The reason for shifting the weighting factor to the results of the load steps was that the focus was on the transient response of the FC and therefore, the main concern was to reproduce the measured load steps with as little error as possible.

Experimental Results
The measurement data used in this work were obtained using a dedicated test stand equipped with a commercial 30 kW fuel cell stack, the associated balance-of-plant components and a specifically designed measurement and control system. As the test stand has also been used for measurements in other publications and to avoid redundancy with other articles, the authors refer the reader to the work in [25,45] for a more detailed presentation of the test stand setup.
For the optimization of the parameters and validation of the simulation results, 6 measurement series were available, divided into three polarization curves and three load steps. Table 1 gives an overview of the measured data series.
For all measurements, the stoichiometry at the cathode was λ cat = 2 for current densities below 0.12 A/cm 2 and λ cat = 1.5 for current densities exceeding 0.22 A/cm 2 . Anode stoichiometry was a product of mixing the recirculated humid hydrogen with dry hydrogen from the hydrogen storage system. The rotational speed of the hydrogen recirculation pump (HRP) was set to a constant value, hence the anode stoichiometry was not constant. It ranged from approx. 6 at low currents to 1.6 at high currents. The anode inlet relative humidity was set to 99% and the pressure was set to a constant value (1700 mbar or 1400 mbar). Ambient air conditions during the measurements can be found in Table 2. Load Step 1700 55 L2 Load Step 1700 55 L3 Load Step 1700 55 A constant air inlet temperature of 40°C could be achieved due to an air cooler downstream of the compressor. The polarization curves were recorded by current controlled stack at different cooling inlet temperatures and anode inlet pressures, as shown in Table 1. For each of the measured polarization curves six measured points in the current range of 0.1-1 A/cm 2 were available, as shown in Figure 2. The series of measurements with the load steps were recorded by voltage-controlled stack. Three different load steps were available: • L1 is recorded with cell potential from 844 mV to 792 mV to 844 mV • L2 is recorded with cell potential from 666 mV to 640 mV to 666 mV • L3 is recorded with cell potential from 834 mV to 630 mV to 834 mV Which corresponds to approximately 0.12 A/cm 2 to 0.22 A/cm 2 to 0.12 A/cm 2 for L1, 0.73 A/cm 2 to 0.85 A/cm 2 to 0.73 A/cm 2 for L2 and 0.12 A/cm 2 to 0.95 A/cm 2 to 0.12 A/cm 2 for L3. Figure 3 shows the measured cell potential and cell current over time.

Hybrid 3D Model Details
The analyzed commercial PEMFC consisted of bipolar plates made of fuel cell grade graphite, carbon paper GDLs with hydrophobic PTFE coating, and an ultra-thin automotivegrade PEM made from SSC-PFSA with reinforcement. The thicknesses of the fuel cell components are listed in Table 3. The model consisted of three different circuits: (1) gas circuit, (2) electrical circuit, and (3) thermal circuit to set the operating temperature. In the direction of the gas flow in the channel, the fuel cell was divided into six slices, which was a compromise between accuracy and computation time. Species mass flows, pressures, and temperatures were applied as described in Section 2.3.

Parameterization
Numerical optimization algorithms were used to find the best values for a selected set of parameters. The following parameters were design variables for the optimization process (1) RECD, (2) activation energy of exchange current density of the CCL, (3) transfer coefficient of the CCL, (4) reference membrane ionic conductivity, (5) activation energy of the ionic conductivity, (6) electrical conductivity, (7) porosity, and (8) tortuosity of the GDL. Table 4 gives an overview of the used material parameters and the design variables are marked with . As the electrochemical parameters of the CCL are mostly unknown, a common practice is to use these parameters for calibration of the electrochemical model, see, e.g., in [18,[46][47][48]. As the temperature influence of the RECD was unknown, the activation energy of the Arrhenius term (see Equation (5)) was selected as a design variable. Although the ionic conductivity could be found in datasheets of the membrane in H-form and for fully hydrated membrane at a reference temperature, the temperature dependence was not available for the given membrane material. Therefore, the activation energy of the temperature dependence and the ionic conductivity were added to the design variables for the optimization process. Due to the lack of information on the water sorption behavior of the used membrane, the model suggested by Springer [5] was applied. The water diffusion coefficient and electro-osmotic drag coefficient were taken from standard ionomer material. For the fluid flow through porous media, the porosity is crucial. As shown by Karpenko-Jereb et al. [22] and d'Adamo et al. [19] the porosity of the GDLs has a significant influence on the species transport, electric conductivity and, therefore, current density. Thus, the porosity was added to the design variables for numerical optimization. Although the tortuosity of the GDLs can be estimated with the porosity as proposed by [49] (Other relationships in the literature can be found in [49][50][51].), it was used as a design variable as well.

Optimization
The results of the optimization process and the comparison between the measured values from the test stand with the simulated results are described in this section. The numerical optimization with the differential-evolution algorithm was started with the bounds shown in Table 5. The differential-evolution strategy was set to 'best1bin'. This strategy chooses two members randomly and their difference is used for the mutation of the best member to construct a trial vector. For more information, see in [40,52]. The initial sampling of the design space was done with the Latin hypercube sequence to maximize the coverage of the available parameter space. The recombination constant (crossover probability) was left to the default number of 0.7. The relative tolerance for convergence was set to 0.01. The global optimization with the differential-evolution algorithm was stopped after 72 h or 1412 function evaluations to limit the computation time. One simulation of a design point takes around 170 s for all six cases. The design point with the lowest integrated difference between simulated and measured value, see Equation 9, was calculated at iteration 957. The set of parameters from the differential-evolution algorithm was used as a starting point for the NMSA to perform a local optimization. The absolute tolerance and function tolerance of the Nelder-Mead Algorithm were set to 0.01. This condition was met after 745 model evaluations. Figure 4 shows the course of the objective value over the iterations of the model. Not every parameter combination leads to a good result. If the simulation time exceeds a certain value, the design point is not feasible and the objective value is set to an artificially high value (1000) not to influence the optimization result. In Figure 4 these values are not shown. The bounds specified in Table 5 were used by the differential-evolution algorithm for the sampling of the design space, whereas the implemented NMSA searched for the minimum without bounds. Therefore, it would be possible to get parameters outside of the specified bounds with the NMSA. The final design variables with the lowest integrated difference are shown in Table 5. Table 5. Bounds of the design variables for global optimization and results from the numerical optimization process. DiffEvo represents the values obtained from the differential-evolution algorithm. These parameters are used as starting vector for the Nelder-Mead Simplex Algorithm to find the local minimum more accurately. The final set of parameters is shown in row NelderMead.

Parameter
i The values obtained from numerical optimization for the design variables agree well with the standard physical values found in the literature, as shown below.

Activation Energies
As in the introduction already mentioned, the measurements showed strong temperature dependence, therefore the activation energies of the Arrhenius terms in the exchange current density (E act,i 0 ,re f ) and the ionic conductivity (E act,σ,ion ) showed higher values than in the literature found. For the activation energy of the exchange current density, a common value is E act,i 0 ,re f = 66 × 10 3 J/mol [22,48] and for the activation energy of the ionic conductivity near E act,σ,ion = 10 × 10 3 J/mol, e.g., E act,σ,ion = 10, 542 J/mol [5] or E act,σ,ion = 9894 J/mol [53] for Nafion ® N117 .

Reference Exchange Current Density and Transfer Coefficient
The exact values for the RECD and transfer coefficient of the CCL were not known. Therefore, these variables were used for the calibration of the model. The value presented, e.g., in [54] with 3 × 10 −9 A/m 2 Pt , for the ORR with Pt-C is specified per unit area of metal, while the model used in this work relates the RECD to the active area of the fuel cell. Since the real surface area of the platinum for the CCL is unknown, the value obtained from the optimization process cannot be directly compared with values from the literature. The optimized transfer coefficient with α = 0.33 lies within the bounds for oxygen reduction reactions (ORR) of hydrogen fuel cells (0.1 to 0.5) as shown in [55].

Ionic Conductivity
The ionic conductivity of the membrane can be found in the datasheet provided by the manufacturer with 14 S/m for fully hydrated membrane at a temperature of 25°C. The obtained value from the optimization process is with 14 S/m at T re f = 55°C lower for the higher reference temperature (The higher the temperature, the higher the ionic conductivity). The reasons for this difference can be that there is an uncertainty if the values from the datasheet were obtained with a fully hydrated membrane activated with 10% H 2 SO 4 at 80°C or the reference water concentration C w,re f = 14 is not the highest possible water concentration of the used membrane. As the water sorption isotherm was not known and the relationship from Springer [5] was used for the simulations, the uncertainty may be in the use of this relationship. However, the obtained value is within physical bounds for SSC-PFSA ionomers [56][57][58].

GDL Electrical Conductivity
The effective electrical conductivity of the GDL σ elec,gdl = 355.8 S/m is higher than given by the datasheet of the manufacturer with 262.5 S/m at 1 MPa compression. The exact value of the stacked cells' compression is unknown, so this could be a reason for the higher value. Another factor is the calculation of the effective value of porous media. The effective electrical conductivity is a function of the bulk electrical conductivity and the porosity of the GDL. Several different relationships can be found in the literature, see in [19]. As proposed by Bruggeman [50] and adapted by Das et al. [59] the relationship of effective and bulk conductivity can be given with With σ bulk the bulk conductivity of the solid and ε the porosity of the material. However, the software used has implemented the following relationship for the effective conductivity : The effective value was specified during the optimization process and the software calculated the bulk electrical conductivity with the above relationship. The bulk electrical conductivity and the solid fraction of the porous media were then used to determine the ohmic losses of the GDLs. The electrical conductivity of commercially available carbon paper GDLs with hydrophobic treatment ranges from 150-350 S/m (effective value in an uncompressed state); thus, the obtained value is comparable.

GDL Porosity and Tortuosity
The porosity with ε GDL = 0.69 and the tortuosity τ GDL = 1.0 obtained from the numerical optimization agree with common values for GDLs in the literature. In addition, the porosity can be estimated with the area weight given by the datasheet from the manufacturer and an estimated compression of approx. 1 MPa with ε GDL,est = 0.7. This value is close to the result of the numerical optimization procedure. The tortuosity of the GDL is mostly unknown, thus several relationships exist in the literature to determine it from the porosity of the material, see in [19]. A value of about τ GDL,expected = 1.38 could be expected from the optimization process by applying the schemes on the determined porosity. However, the obtained value of τ GDL = 1.0 is low. The diffusion length of the GDL and the thickness are not the same in reality. Moreover, τ GDL = 1.0 is equivalent to free diffusion of gas species through the GDL as the binary and actual diffusion coefficients, see Equation (3) have the same value.
With the parameter set from the numerical optimization the following results can be obtained. The simulations showed very good agreement with the measurements, as can be seen in Figure 5. Although the weighting factors favored the load levels, the results of the simulated polarization curves differed only slightly from the measured data with average relative errors of ∆P1 = 1.01 %, ∆P2 = 1.35 %, and ∆P3 = 0.54 %. The load step in the high current region, L2, agreed very well, whereas L1 in the low current region showed more deviation to the comparison data. However, the load step over the entire current range, L3, agreed well with the measurements. In addition, the current after the load step was represented well by the model. The small current steps before and after the main load steps came from the measurement from the test stand and were necessary to ensure the stable operation of the stack.  Figure 6. Negative values indicate mass flows at outlets. It can be observed that the water vapor mass flow at the anode outlet in Figure 6a increased after the load step and the operation of the fuel cell at 0.95 A/cm 2 . The relative humidity also increased and remained at a constant value until the current dropped back to 0.12 A/cm 2 . The relative humidity at the cathode outlet approached a maximum value of 1.5 after the load step to 0.95 A/cm 2 . The water vapor mass flows at the anode and the cathode outlet increased steadily. Due to water diffusion through the membrane, the water vapor concentration at the anode also increased.  Figure 7 shows the membrane water content C w for all three load steps over time. C w for load step L1 is lower than for L2 and L3. The current densities L2 and L3 after the load step were much higher than for L1. Therefore more water was produced by the chemical reaction in the CCL. This water increased the relative humidity in the cathode and was sorbed by the membrane. The maximum membrane water content was given with 17.5.

Figure 7.
Membrane water content C w for load steps L1, L2 and L3. C w depends on relative humidities at interface to anode GDL and CCL. Max. membrane water content was C w,max = 17.5

Conclusions
The parameterization of mechanistic models for PEMFCs is challenging. Developing more sophisticated models with, e.g., anisotropic material properties, agglomerate models for CL or gas crossover in the membrane further exacerbates the problem. The measurement methods to determine the necessary parameters are complex, expensive, and specialized equipment is often needed, e.g., the porosity of the GDL can be measured via 3D XCT methods, see [19]. This work demonstrated the efficient parameterization of a single-phase, isothermal PEMFC model with measurements on a stack. The measurements showed a strong temperature dependence and the operating temperature is relatively low with 52-55°C (cooling inlet). Therefore, special attention was paid to the activation energies of the Arrhenius terms of the RECD and the ionic conductivity of the membrane.
The following parameters were selected as design variables for numerical optimization: (1) RECD of CCL, (2) activation energy of the exchange current density, (3) transfer coefficient of the CCL, (4) reference membrane ionic conductivity, (5) activation energy of the ionic conductivity, (6) electrical conductivity, (7) porosity, and (8) tortuosity of the GDL. The numerical optimization was conducted with the Python package scipy.optimize. The differential-evolution algorithm for global optimization and the Nelder-Mead algorithm for local optimization were used. It was shown that both algorithms were well suited for this type of engineering problem and fast convergence of the objective value was given. With these optimized parameters, the simulated and measured polarization curves and load steps agreed very well. Although numerical optimization had found good values for the selected design variables, it does not necessarily replace advanced measurement methods to determine the unknown parameters in advance. Further activities will be necessary here in the future.
However, the measurement of the load steps and polarization curves in conjunction with the in-and outlet mass flows, pressures, humidities and temperatures, and cooling temperatures were sufficient to parameterize the fuel cell model used. Due to the lack of measured material data, the optimized parameters were compared with values from the literature. They agreed well with the typical physical values of the eight design variables, although the exact values remain unknown. The obtained material parameters will be used in the future to parameterize a more detailed 3D CFD model to analyze the transient behavior of PEMFCs at low operating temperatures with high spatial and temporal resolution.

Acknowledgments:
The computational results presented have been achieved, in part, using the Vienna Scientific Cluster (VSC). The authors acknowledge TU Wien Bibliothek for financial support through its Open Access Funding Programme.

Conflicts of Interest:
The authors declare no conflicts of interest.