Data-Driven Reduced-Order Modeling of Convective Heat Transfer in Porous Media

: This work presents a data-driven Reduced-Order Model (ROM) for parametric convective heat transfer problems in porous media. The intrusive Proper Orthogonal Decomposition aided Reduced-Basis (POD-RB) technique is employed to reduce the porous medium formulation of the incompressible Reynolds-Averaged Navier–Stokes (RANS) equations coupled with heat transfer. Instead of resolving the exact ﬂow conﬁguration with high ﬁdelity, the porous medium formulation solves a homogenized ﬂow in which the ﬂuid-structure interactions are captured via volumetric ﬂow resistances with nonlinear, semi-empirical friction correlations. A supremizer approach is implemented for the stabilization of the reduced ﬂuid dynamics equations. The reduced nonlinear ﬂow resistances are treated using the Discrete Empirical Interpolation Method (DEIM), while the turbulent eddy viscosity and diffusivity are approximated by adopting a Radial Basis Function (RBF) interpolation-based approach. The proposed method is tested using a 2D numerical model of the Molten Salt Fast Reactor (MSFR), which involves the simulation of both clean and porous medium regions in the same domain. For the steady-state example, ﬁve model parameters are considered to be uncertain: the magnitude of the pumping force, the external coolant temperature, the heat transfer coefﬁcient, the thermal expansion coefﬁcient, and the Prandtl number. For transient scenarios, on the other hand, the coastdown-time of the pump is the only uncertain parameter. The results indicate that the POD-RB-ROMs are suitable for the reduction of similar problems. The relative L 2 errors are below 3.34% for every ﬁeld of interest for all cases analyzed, while the speedup factors vary between 54 (transient) and 40,000 (steady-state).


Introduction
The simulation of fluid flows is essential in many engineering fields. For this, highfidelity models using Computational Fluid Dynamics (CFD) are commonly employed to numerically solve the Navier-Stokes equations on two-or three-dimensional domains. These CFD models are often computationally expensive; therefore, they are not adapted for multi-query applications, such as uncertainty quantification or design optimization, where numerous solutions are needed for various model/input parameter values. Model-Order Reduction (MOR) techniques can lower the computational burden of such applications by the generation of Reduced-Order Models (ROMs), which are considerably faster to solve. Thermal-hydraulic system codes serve as good examples for ROMs since the Navier-Stokes equations are simplified to one-or zero-dimensional equations by solving for surface or volume averaged values of the variables [1][2][3]. These reduced models are often referred to as physics-based ROMs and are commonly used for the analysis of industrial systems [4][5][6]. A disadvantage of these codes is that they often struggle with the modeling of multi-dimensional effects in larger fluid volumes (flow detachments, recirculation, etc.). zones. Lastly, Section 4 discusses the results, while Section 5 draws our conclusions and proposes further improvements for the method in this paper.

Materials and Methods
In this work, a Proper Orthogonal Decomposition aided Reduced-Basis (POD-RB) technique is developed for the parametric Model Order Reduction (MOR) of hydraulic systems where the structural elements are homogenized and treated as porous media. The flow chart of the MOR process is presented in Figure 1. It starts with the identification of the governing laws for the fluid flow on a specific domain with appropriate boundary conditions.
Then, we enter the training phase, which requires the generation of a Full-Order Model (FOM) by the discretization of the previously selected governing equations. Following this, the FOM is exercised to collect information about the parameter-and time-dependent solution manifold. The acquired information is then synthesized in reduced operators. This concludes the training phase, which is computationally expensive since it requires several solutions of the FOM.
The process continues with the online phase, which consists of the rapid assembly of the Reduced-Order Models (ROMs) using the precomputed reduced operators, the solution of the ROM, and the evaluation of chosen Quantities of Interest (QoIs). The steps of the MOR process are discussed in detail in the following subsections.

Identifying the Problem (Governing laws and domain)
Creating

Governing Equations and the Full-Order Model
In many cases, the high-fidelity flow configuration of systems with structural elements cannot be resolved explicitly due to a prohibitively high number of cells needed, the difficulty to generate a mesh for the complex geometry, or limited computational resources, etc. Therefore, these complex structural elements are often homogenized into porous medium zones, and the interactions between the fluid and the homogenized solid structures are modeled using several semi-empirical friction correlation functions. The derivation of the momentum and mass conservation equations (Navier-Stokes equations) in the porous media was discussed in [7], and only the final form is presented here.
The pressure jumps at the porous medium interfaces are not modeled in this work; we assumed that these jumps are negligible compared to the pressure losses within the homogenized structures. This simplification should not significantly impact the flow field since the thermophysical properties of the flow are pressure independent. Nevertheless, if needed, concentrated pressure drops at the inlet and outlet could, in principle, be simulated by adding additional inlet and outlet regions with suitable momentum sinks. We denote the computational domain by Ω and the corresponding boundary by Γ. Considering this notation, the porous medium conservation equations for Newtonian fluids can be expressed as: where u D = u D (r, t) is the Reynolds-averaged Darcy velocity vector field, p = p(r, t) is the Reynolds-averaged corrected pressure field, T = T(r, t) is the Reynolds-averaged temperature field, η is the molecular dynamic viscosity, η t is the turbulent eddy viscosity, ρ is the density, and β th is the thermal expansion coefficient. In this work, ρ, η, and β th are considered to be constants. Furthermore, g denotes the gravitational acceleration, γ = γ(r) is the volume fraction occupied by the fluid, while F p = F p (r, t) and F f r = F f r (r, t) are volumetric linear momentum sources and sinks resulting from the interaction between the fluid flow and the homogenized structural elements. The Darcy velocity can be expressed as where u(r, t) is the real velocity vector field. In most of the CFD solvers generated using the OpenFOAM © framework [35], the utilization of the corrected pressure (p) is preferred over the total pressure (p t ). The following relationship between p and p t can be derived using Green's identity and the fact that the density, thermal expansion coefficient, and the gravitational acceleration are considered to be constant fields [36]: where T ref denotes a reference temperature. Note that, in the case of clear fluid regions (γ = 1), Equations (1) and (2) turn into the incompressible Reynolds-Averaged Navier-Stokes (RANS) equations [37] with Boussinesq eddy viscosity [38] and buoyancy approximations [39]. For the closed-loop problems in this work, the boundary conditions on the wall segments (Γ wall ) for the fluid dynamics fields can be expressed as: In case of transient simulations, we assume that the initial velocity (u D,0 ) and pressure (p 0 ) fields are given and satisfy the boundary conditions. Due to the fact that no Dirichlet boundary conditions are specified for the incompressible fluid flow problem, and that Equation (2) contains only ∇p, the corrected pressure in the system cannot be uniquely determined. The allowable pressure solutions are shifted by a constant pressure field.
For this reason, in the closed-loop problems in this work, we impose the pressure at a specific location (r 0 ) to set the correct operating point: p(r 0 , t) = p 0 . Terms F f r and F p in the momentum equation describe the linear momentum sources and sinks resulting from the interactions between the fluid and the homogenized structures in the porous medium. We assume that the momentum source, F p , corresponds to a homogenized pump, which can be modeled by a vector field with region-wise constant values in the N pump pump regions (Ω pump ) and zero everywhere else: The flow resistance, on the other hand, cannot be considered independent of the fluid velocity. Based on [7,40], the following empirical correlation function is used to approximate the obstructing force in direction i = {x, y, z}: where D h is the the hydraulic diameter of the structure and f denotes the friction factor, which can be computed using the following expression: where the A f and B f parameters can be determined considering the flow regime or the geometry of the problem. Lastly, we note that η t can be computed using different eddy viscosity models. In this work, a k − model is used, which requires the solution of two additional equations for the turbulent kinetic energy (k) and its dissipation rate ( ). For more details on this turbulence model in clean fluids and in porous medium flows, see [7,41], respectively. The porous medium momentum and mass conservation equations are supplemented with a porous medium enthalpy equation, which is solved for the temperature field in the system: where c p and k l are the specific heat and thermal conductivity of the fluid. We assume that c p and k l are constant parameters. The second term on the right-hand side of Equation (9) describes a heat-exchange process between the fluid and the structural elements. This is characterized by h, the heat transfer coefficient, T ext , the temperature of the structural element, and A V , the volumetric fluid contact area of the structural element. Term Q, on the other hand, describes a volumetric heat source in the system, which can result from chemical or nuclear reactions, to name two examples. Moreover, the additional mixing effects from the turbulence are taken into account using a turbulent diffusion term with α t = η t Pr t diffusivity, where Pr t denotes the turbulent Prandtl number. For simplicity, we assumed that the closed-loop systems in this work are equipped with proper thermal insulation; therefore, zero gradient boundary conditions can be employed for the temperature field: ∇T · n = 0, r ∈ Γ.
For transient simulations, we expect that the initial temperature field T 0 is known and satisfies the boundary conditions.
Lastly, it is important to note that the governing equations are parametric in the sense that several model parameters, such as the material properties or parameters of the homogenized structures, may be uncertain. We denote the vector containing these uncertain model parameters by µ = [µ 1 , ..., µ N µ ] T .
Once the governing equations are identified, the process continues with the offline phase, which starts with the generation of a Full-Order Model (FOM). In this work, this is handled by GeN-Foam [7] (GeN-Foam repository: https://gitlab.com/foam-for-nuclear/ GeN-Foam, accessed on 12 August 2020), an OpenFOAM © [35,42]-based multiphysics solver. The spatial discretization in OpenFOAM © employs the cell-centered Finite Volume Method (FVM) [42,43], which is popular in many CFD applications. This process begins by splitting the spatial domain Ω into N non-overlapping cells, Ω i , (i = 1, ..., N), with volumes of V i . Then, the governing equations are integrated over each cell with the assumption that the solution is constant within the cells and on the cell boundaries. Given that the equations are linearized, the integration results in a linear algebraic equation system, which can be solved using an iterative method. For more details on how specific differential operators are discretized, we refer the reader to [42]. The discretization in time is handled using an implicit Euler scheme in this work, while the nonlinearities are resolved by fixed-point iteration.

Training Reduced-Order Models
The training of the ROMs continues with the collection of information about the parameter-and time-dependent solution manifold. This is carried out using the method of snapshots [11], in which the FOM is exercised to obtain N s snapshots of the system at different time instances and/or model parameter values. In this work, we assume that N µ different parameter samples are used and snapshots are saved at N τ different time instances, yielding a total number of N s = N µ × N τ snapshots for a parametric transient and N s = N µ for parametric steady-state runs.
For FOMs built using FVM, a snapshot is a vector of length N or 3N (assuming three spatial dimensions) for scalar and vector fields, respectively. The snapshots for field X = {u D , p, T, F f r,z , η t , α t } are stored in corresponding snapshot matrices denoted by R X : where we see that, in addition to the snapshots of the solution fields, we capture the eddy viscosity (η t ), eddy diffusivity (α t ), and the flow resistance force in porous medium zone z = 1, ..., N z as well. The importance of saving these snapshots is discussed later in detail. Next, we introduce the reduced-basis approximations for the fields of interest. We assumed that the time-and parameter-dependent solution vectors live in a considerably smaller subspace of R N . In other words, we assume that the fields of interest can be approximated as the linear combination of only a few, space-dependent global basis vectors (ψ X ). This also means that the summation coefficients (c X = c X (t, µ)) in the linear combination are responsible for capturing the time-and parameter-dependence of the approximate solution. For the Darcy velocity and corrected pressure fields, this can be expressed as follows: Instead of the real velocity, we approximate the Darcy velocity (u D ). The reason for this is that u D is smooth in space and, therefore, can be better approximated using global basis functions. Similarly to the velocity and pressure, the temperature field can be approximated as: Moreover, for simulations involving turbulence, the turbulent eddy viscosity (η t ) and the turbulent diffusivity (α t ) are also approximated as follows: Lastly, a reduced-basis approximate is constructed for the flow resistance in porous medium zone z as: In this work, the spatial basis vectors of the fields of interest are obtained via Proper Orthogonal Decomposition (POD) [10,13,44] applied to the snapshots. This involves the discovery of a reduced set of orthonormal basis functions that span a subspace of fixed rank that is closest to the snapshots in a global L 2 sense. For field X, this optimization problem can be expressed as follows: where a, b Ω denotes the inner product between fields a and b. For FVM, this equivalent to following matrix product a T M Ω b where M Ω is a mass matrix containing the volumes of the mesh cells on the diagonal. In practice, the basis functions of the reduced subspaces can be obtained using a constrained Singular Value Decomposition (SVD) [45] of the snapshot matrix: Once the SVD of R X is computed, the first r X columns of the left singular vector matrix Ψ X can be used as the basis vectors for the reduced-basis approximation. The selected basis vectors are often referred to as POD modes, and the two expressions are used interchangeably in this work. Moreover, r X is the rank of the subspace and can take a maximum value of N s . In the examples presented in Section 3, an energy truncation strategy is used to determine the actual value for r X . This utilizes the singular values (σ X i , i = 1, ..., N s ) on the diagonal of Σ X : where τ is a truncation threshold. In general, as the truncation value decreases (i.e., as 1 − τ → 1), more basis functions are used to build the reduced subspaces. The basis functions for X = {u D , p, T, F f r,z , η t , α t } are constructed using this approach.
Once the basis functions are computed, the reduced-basis approximations in Equations (12)-(15) are plugged back into Equations (1), (2) and (9). Following this, the equations are projected using the basis functions of the corresponding solution variable (Galerkin projection). This results in the following reduced equation system for the linear momentum and mass conservation equations: where the elements of the low order vectors, matrices, and higher-order tensors are computed as: where δ z (r) is a selection function, which returns 1 if r is in porous medium zone z and 0 otherwise. In this scenario, the expansion coefficients of the temperature (c T ) are considered to be fixed. To be able to solve system (19) and (20) for the expansion coefficients of the velocity and pressure (c u D , c p ), however, three additional issues need to be addressed: 1.
In this form, the two-equation ROM is not stable since the Ladyzhenskaya-Babuska-Brezzi (LBB) condition [46] is not necessarily satisfied. For this reason, an approximate supremizer stabilization [14,18] is used. This entails the generation of supremizer snapshots using the pressure snapshots by solving the following problem: The supremizer snapshots are collected in a snapshot matrix (R s ) and the corresponding Ψ s basis functions are generated using the constrained SVD. Then, the velocity reduced subspace is augmented by these basis functions. In practice, this means that the supremizer basis vectors are appended to the velocity basis matrix: The coefficients of the eddy viscosity and diffusivity are still unknown. The computation of these coefficients from reduced forms of eddy-viscosity models is not widely utilized by the fluid dynamics ROM community; therefore, a non-intrusive, data-based approach-Radial Basis Function (RBF) interpolation-was chosen in this work. The effectiveness of this method has been demonstrated in [19][20][21]34]. This approach interpolates coefficients c η t and c α t within the parameter space and time using the snapshots as anchor nodes. 3.
The coefficients of the flow resistances are still unknown. In this work, the Discrete Empirical Interpolation Method (DEIM) [47] is used for the determination of these models based on the coefficients of the Darcy velocity. This can be expressed as: where P T z is a matrix that selects r F f r,z rows when applied to a vector. In this work, P T z was constructed using the algorithm published in [47]. This method ensures that the evaluation of the flow resistance at reduced-order level is computationally efficient.
As a next step, the porous medium enthalpy equation is reduced with the same process. The reduced-basis approximates are plugged into Equation (9), and the resulting expression is projected using the basis functions of the temperature: where the elements of the reduced operators and source terms can be expressed as:

Evaluating Reduced-Order Models
Once the training phase of the ROM is concluded, we enter the online/evaluation phase. This phase starts with the assembly of the ROM using a new parameter sample (µ * ) where the quantity of interest needs to be evaluated. This is computationally cheap if the reduced-operators are affine in the model parameters and time, which allow us to express the reduced-order operators as a sum of products of space-dependent operators (A i ) and solution-, parameter-, and time-dependent scalar functions ( f i (c X , t; µ)) as: In this scenario, the assembly of the reduced-order model simplifies into the summation and scalar multiplication of small, precomputed matrices and vectors from the training phase. By looking at Equations (19), (20) and (24), we conclude that the ROMs are affine in the following parameters: the density (ρ), the magnitude of the pumping power in the porous medium zones (|F p,z |), the thermal expansion coefficient (β th ), the molecular dynamic viscosity (η), the specific heat (c p ), the thermal conductivity (k l ), the heat transfer coefficient (h z ), and the external coolant temperature (T ext ) in the heat exchangers. Furthermore, through the utilization of DEIM, the problem becomes affine in the parameters of the flow resistance function as well.
Since the ranks of the selected subspaces are typically much lower than the size of the original, full-order system (r X N), the assembly of the ROM is computationally cheap. This process is followed by the solution of the ROM for the expansion coefficients. In this work, an implicit Euler discretization is used for time-dependent problems, while the nonlinearities are resolved using a nested fixed-point iteration between the fluid ROM and enthalpy ROM. Once the expansion coefficients are obtained, the reduced-basis approximates can be reconstructed and the quantities of interest (point values, averages, etc.) can be evaluated at the new parameter sample (QoI(µ * )).

Results
We demonstrate the applicability of the developed ROM to parametric steady-state and transient problems in the Molten Salt Fast Reactor (MSFR) [33]. For this, a 2D axisymmetric model is employed. The MSFR is a liquid-fuel nuclear reactor in which the fissile material is dissolved in a molten salt that is circulated in a loop. Heat (3 GW th ) is generated by the nuclear fission reactions and is predominantly deposited in the liquid fuel. Figure 2 shows the geometry and mesh used for the simulations together with the distribution of the power deposition. The mesh was adopted from [34] and contains 16,140 cells. For every scenario in this work, this distribution is considered to be fixed. The red and blue regions in the mesh denote porous medium zones, while clean fluid is shown in green. The molten salt fluid enters the core cavity through the core inlet, warms up, and leaves through the core outlet. The molten salt then cools off in the heat exchanger. This example has been specifically chosen to showcase how a system with both clean fluid and prorous medium zones can be treated with the developed ROM.
Even though we use a nuclear reactor as an example, the developed ROMs can be used for chemical reactors and other industrial loops involving heat transfer. Both steadystate and transient examples assume a turbulent flow regime, with approximate Reynolds numbers in the [6.72 × 10 5 , 1.23 × 10 6 ] interval within the core cavity. The default values of the thermophysical parameters of the fluid were determined using the help of [48] and are summarized in Table 1. As already mentioned, the pump and the heat exchanger can be treated as porous medium zones. In this work, the pump is treated as a volumetric momentum source with the assumptions that the fluid entirely occupies the volume (γ = 1), and the force exerted by the pump is homogenized in the region with a 100 kN m 3 default value. This is justified by the lack of data on the friction factors and structure of the pump in the MSFR design. The heat exchanger, on the other hand, has been treated as a porous medium with default parameters listed in Table 2.
These parameters are considered to be approximate values, since no specific design has been chosen for the MSFR yet. The parameters describing the flow resistance have been chosen using the Blasius correlation in [49] (assuming a tubular heat exchanger) with slightly increasing A f to obtain a realistic pressure drop in the porous medium zone. External temperature T ext 850 K The first numerical example considers a steady-state scenario of the MSFR with several uncertain parameters. As the second case, a transient scenario is analyzed. The transient starts from a steady-state, which is obtained using the default parameters from Tables 1 and 2. At t = 2 s, a flow coastdown is simulated by linearly ramping down the pumping force from 100 kN m 3 to 20 kN m 3 in ∆t ramp . Figure 3 showcases this transience with ∆t ramp = 8 s.
The problem at hand was discretized in space using corrected linear gradient estimators for the diffusion terms (second order) and a vanLeer [50] scheme for the advection term in the fluid-dynamics equations (also second order). The advection term in the enthalpy equation was discretized using a first order upwind scheme. In terms of time discretization of transient scenarios, both the FOM and the ROM applied a backward Euler scheme.

Reduced-Order Models for Steady-state Problems
First, the applicability of the developed ROM to steady-state scenarios was assessed. For this, five parameters were assumed to be uncertain. The list of these parameters together with the assumed probability distributions is presented in Table 3. In this context, U (a, b) is a uniform distribution in the [a,b] interval. The generation of the ROM for the parametric steady-state simulations starts with collecting snapshots of the fields of interest. Altogether, N s = 40 snapshots were generated by sampling parameter vector µ = [|F p |, T ext , h, β th , Pr] T using Latin Hypercube Sampling (LHS), executing the FOM using the generated samples, and saving the given fields into snapshot matrices. The number of training samples was determined considering the number of uncertain parameters, the available time/computational resources, and resolution constrains for the RBF interpolation of the turbulent eddy viscosity/diffusivity coefficients.
Then, the supremizer snapshots were obtained by solving Equation (22). Following this, the basis functions of the reduced subspaces were extracted from the snapshots using a constrained SVD. The first 3 POD modes of the velocity, pressure and temperature are presented in Appendix A. The decay in the squares of the singular values of the snapshot matrices is presented in Figure 4. Note, the values were normalized by the square sum of the singular values. Furthermore, since the energy truncation limit (τ) was employed in this example for the determination of the number of POD modes used for the generation of the ROMs (see Equation (18)), the truncated energy per POD mode is showcased in Figure 4 as well.
Overall, we see a rapid decay in the squared singular values, indicating that the solution exhibits a low-dimensional behavior.
Then, numerous ROMs were prepared with different numbers of basis functions. The number of basis functions for the fields of interest was determined using Equation (18), with two exceptions. The number of modes for the pressure and supremizer fields was determined in a slightly different manner since the inclusion of too many supremizer modes may decrease the accuracy of the ROM. In this work, we used the same number of pressure and supremizer modes. Once the number of POD modes for every other field was chosen using the truncation limit, we performed a one-dimensional optimization to select the number of pressure/supremizer modes, which yielded the highest average accuracy in velocity over the test set.
Furthermore, due to the limited number of snapshots and the relatively high-dimensional parameter space, the RBF interpolation for coefficients of the eddy viscosity and diffusivity is restricted to the pumping-force alone, simplifying it to a one-dimensional problem. This can be justified by the fact that the other parameters influence the fluid dynamics equations through the buoyancy effect, which is significantly lower than the effect of changing the pumping force itself.
As a next step, a validation set was assembled by drawing N val = 20 new parameter samples using LHS. The size of the test set was determined to be half of that of the training set due to the FOM being relatively expensive to solve, limiting the available data. These samples were used to test (validate) the ROMs with respect to the FOM. The parameter vectors in the validation set were plugged back into the FOM, and the different ROMs and the corresponding FOM solution fields were compared to the reconstructed reduced-basis approximates from the ROMs. We chose the relative L 2 error for an indicator: where X = {u D , p, T, η t , α t }. The average and maximum errors over the validation set as functions of the truncation parameter τ are presented in Figure 5. We see a considerable decrease in the error as we decrease the truncation parameter. The error in the eddy viscosity went to saturation. This was caused by the fact that the error in the RBF interpolation became more dominant compared to the approximation error beyond τ = 10 −4 . The behavior of α t was exactly the same as η t since the two fields (in this work) differ only by a constant factor (turbulent Prantl number, Pr t ). The ROM built with τ = 10 −8 was selected for further investigation. The number of modes used for the fields of interest together with the error statistics over the validation set are presented in Table 4. Moreover, the selected ROM utilized seven modes for the flow resistance in the heat exchanger during the DEIM process. It is visible that, with the exception of η t /α t , the maximum relative L 2 error was below 0.53%, which is considered to be adequate. Since the relative L 2 error is an integrated quantity, we present the distribution of the local errors of the velocity and temperature fields at the test sample, which exhibited the highest errors in terms of the velocity. Figure 6 shows the FOM solution (to the left) and the absolute error between the FOM and the ROM (to the right). It is visible that the maximum error was attained at the outlet of the heat exchanger within the corner element in the mesh. At this point, the relative error was close to 5% locally. In other parts of the reactor, we see considerably lower relative errors. Note that the error appears to be especially low in the heat exchanger meaning that the DEIM process seems to be efficient for the approximation of the flow resistance.   We see that the maximum error was attained at the bottom of the core cavity. At this point, the relative error slightly exceeded 1% locally, while, at other parts of the reactor, the error remained significantly lower.
Lastly, we evaluated the computational time savings of using a ROM (built with τ = 10 −8 ) instead of the FOM. Obtaining the FOM solutions for a parameter combination took approximately 4000 s on average, while the same was 0.1 s using the ROM. This means that an approximate single-run speedup of 40,000 was achieved.

Reduced-Order Models for Transient Problems
The transient scenario used in this work simulates a pump coastdown from 100% to 20% of the nominal power. The transient is visualized in Figure 3. We consider the time needed to reach the new power level an uncertain parameter with a uniform distribution: ∆t ramp ∼ U (3, 13) s. The generation of the reduced-order model starts with the collection of snapshots. Altogether, five transients are solved evenly spaced in coastdown times. The number of training samples was chosen to ensure an adequate accuracy in the RBF interpolation of the eddy viscosity/diffusivity coefficients while still considering the relatively expensive solutions of the FOM.
The training scenarios are presented in Figure 8. It is visible that there was an oscillation in the maximum temperature in the system, which depended on the dynamics of the velocity field. We took snapshots of the system every 0.25 s for each scenario meaning that, altogether, 805 snapshots were used to build the reduced subspaces. The decay in the normalized squared singular values and the truncated energy are presented in Figure  9. We see that the singular values decayed much slower compared to the steady-state results. This implies that the time-dependence in the problem introduced considerable spatial variation in the fields of interest. As a next step, the basis functions of the fields of interest were built using POD. The first 3 POD modes of the velocity, pressure and temperature are presented in Appendix B. The number of basis functions were determined using the strategy already discussed in Section 3.1, which involves the definition of an energy truncation limit (τ). Several ROMs were then prepared with different values of τ. These ROMs were tested on 10 additional pump coastdown times that were randomly sampled from the [3,13] s interval (test set). In must be noted that, in this scenario, the ratio of test samples to training samples is higher because the expected variation of the solutions in time and parameter space is larger compared to the steady-state problem.
The relative L 2 error defined in Equation (30) was evaluated for each solution field over time. The maximum and average errors over the test set and time as function of the energy truncation parameter are presented in Figure 10. We see a decrease in the average error in the 10 −2 − 10 −7 interval. Beyond this, the errors slightly increased, which can be explained by the buildup of the errors discussed later in Section 4.
The number of modes together with the average and maximum errors for τ = 10 −7 are presented in Table 5. We see that over all the solution fields, times, and test sets, the maximum relative L 2 error was 3.3%, which is adequate for engineering applications. Furthermore, we note that the ROM was capable of the reconstruction of the temperature field with a maximum L 2 error slightly below 0.25%, which is also considered to be accurate. The number of basis functions used to build this ROM is considerably higher than those presented in the steady-case example. This is the direct consequence of the large spatial variation of the solution fields introduced with the time-dependence. The time-dependent error curves of the field variables over the test parameter set are presented in Figure 11 for τ = 10 −7 . It is visible that, for the velocity, pressure and temperature, the highest relative errors appeared in the new steady-state in many cases. The possible reasons behind these effects are further discussed in Section 4.
Lastly, the savings in computation time were evaluated. The execution of the FOM took approximately 16,000 s on a single core of a processor, while the same was around 310 s for the ROM built using τ = 10 −7 . This yielded a moderate speedup factor of 51.6. The relatively slow ROM was the result of one of the shortcomings of the RBF interpolation, namely that the 41 coefficients (for both α t and η t ) need to be interpolated using 805 anchor points at every single time step. Other operations needed for the assembly and solution of the ROM are of considerably lower operation count. Assuming that using 21 POD modes are sufficient for the eddy viscosity/diffusivity (corresponds to τ = 10 −6 ), the execution time of the ROM decreased to 170 s on average, which resulted in a speedup factor of 94.1.

Discussion
Based on the results presented in Section 3, several observations were made about the performance of POD-RB-ROMs for porous medium flows.
First, the supremizer stabilization method for the ROM was not optimal for the handling of turbulent buoyant flows described using Equations (2) with the corrected pressure. This may be the reason behind the ROMs yielding higher relative errors at the new equilibrium levels in the transient scenarios where the buoyancy forces are considerable compared to the momentum source of the pump. Originally, the supremizer stabilization was derived for steady-state clean laminar fluid flows [14]; therefore, the presence of spacedependent diffusivities, porosity factors, and buoyancy-altered corrected pressure terms may require the modification of the original concept. A possible approach could be the inclusion of these aspects in the supremizer equation in Equation (22).
Second, it is visible that the errors for transient problems may slightly increase if too many POD modes are used (beyond τ = 10 −7 ). This can be the result of the utilization of first-order upwind schemes for the treatment of the advection term in the enthalpy equation, which leads to numerical inconsistency between the FOM and the ROM. The reason behind this is that the discretization of the ∇ · (ψ u D i ψ T j ) terms depends on ψ u D i itself, which is not necessarily the same as the discretization used for ∇ · (u D T). This effect can be further amplified by the inclusion of supremizer modes, which do not have a physical meaning. Third, the DEIM-based method for the handling of the interactions between the fluid and the homogenized structure showed promising results. Even though the relative errors in the flow resistances were not evaluated separately, the low errors in the velocity and pressure ROMs indicate that the flow resistances were treated properly.
Fourth, we observed that the adopted RBF interpolation technique had shortcomings for both steady-state and transient problems. In the case of parametric steady-state scenarios, the parameter space may have to be shrunk to obtain adequate accuracy for the interpolation since, in higher parameter spaces, with a limited number of anchor points, the method yielded results with low accuracy. However, neglecting parameter dimensions in the interpolation introduced additional errors as well.
In the case of transient scenarios, the RBF interpolation may considerably slow down the execution of the ROM due to the fact that some transient scenarios may require a large number of POD modes for the eddy viscosity/diffusivity. Coupled with a large number of anchor points in the interpolation, this may yield relatively slow ROMs. A possible improvement on this part may be the utilization of the sparsity in the problem in a sense that the interpolation at a given parameter sample should only evaluate contributions from anchor points that are in its vicinity.

Conclusions
A data-driven Reduced-Order Model (ROM) was developed for the incompressible porous medium RANS equations coupled with the porous medium enthalpy equation. The ROMs were generated using an intrusive Proper Orthogonal Decomposition aided Reduced-Basis (POD-RB) technique coupled with the method of snapshots. The process involved the exploration of the parameter-and time-dependent solution manifolds by solving a higher-fidelity model (or Full-Order Model, FOM) and saving the solution fields as snapshots. The information in these snapshots was then compressed using POD into global spatial basis functions, which were then used to approximate the higher-fidelity solution and to project the equations onto reduced subspaces.
This resulted in a considerable reduction in the spatial degrees of freedom, and, consequently, in the execution time as well. The developed ROM utilized a supremizer stabilization for the reduced porous medium fluid dynamics equations and a Radial Basis Function (RBF) interpolation for the handling of the turbulent viscosity/diffusivity at reduced-order level. The flow resistance exerted by the homogenized structural elements in the flow was treated using the Discrete Empirical Interpolation Method (DEIM).
To the best of the authors' knowledge, the developed ROM is a novel contribution to the fluid dynamics ROM literature. Moreover, this method can be considered as an alternative for approaches where fluid dynamics POD-RB-ROMs are used for large volumes of clean fluids and thermal-hydraulics system codes for other elements in the loop [32].
The ROM has been tested using a 2D axisymmetric model of the Molten Salt Fast Reactor (MSFR), which contains clean fluid and porous medium zones in the same thermalhydraulic loop. The results indicate that the developed ROMs can be used to accurately emulate parametric steady-state and transient problems. The relative L 2 errors of the POD-RB-ROMs were always below 3.3%, which is considered adequate, especially for engineering applications. While the achieved speedup is on the order of 10,000 for the steady-state scenarios, it decreases to 51-94 for transient cases. This can be explained by the shortcomings of the RBF interpolation.
Overall, from an engineering perspective, the developed ROMs yielded accurate solutions. Compared to the accuracy requirements for steady-state simulations with thermal-hydraulics system codes in the nuclear industry, the reduced-order model errors were well below the recommended limits [51].