A Three-Dimensional Homogenization Approach for Effective Heat Transport in Thin Porous Media

Heat transport through a porous medium depends on the local pore geometry and on the heat conductivities of the solid and the saturating fluid. Through upscaling using formal homogenization, the local pore geometry can be accounted for to derive effective heat conductivities to be used at the Darcy scale. We here consider thin porous media, where not only the local pore geometry plays a role for determining the effective heat conductivity, but also the boundary conditions applied at the top and the bottom of the porous medium. Assuming scale separation and using two-scale asymptotic expansions, we derive cell problems determining the effective heat conductivity, which incorporates also the effect of the boundary conditions. Through solving the cell problems, we show how the local grain shape, and in particular its surface area at the top and bottom boundary, affects the effective heat conductivity through the thin porous medium.


Introduction
Heat conduction in porous media is a relevant process in applications ranging from geothermal engineering to various technical applications. Especially in the latter field, many porous components have a thin shape (see, e.g., Belgacem et al. 2016, Michaud 2016, where filters, fuel cells and membranes count among typical examples. A thin porous layer is often part of composite materials (see, e.g., Asbik et al. 2006). For a porous medium where grains and the saturating fluid are under local thermal equilibrium, the effective heat conductivity of the porous medium characterizes the medium's ability to transport heat via conduction. A better understanding of the heat transport in porous media, in particular finding the medium's effective heat conductivity, can help to not only predict the heat transport in a certain setup, but also to design porous materials according to the needs of industrial applications. In this context, the detailed properties of the considered porous medium are necessary to investigate (see, e.g., Ranut and Nobile 2014).

3
When investigating heat transport in porous media, completely pore-scale resolved models (see, e.g., Koch et al. 2021) are often not feasible for large domains due to their computational complexity. Because of that, we will determine the effective heat conductivity of the porous medium at the Darcy scale. While simple approaches like a porosityweighted averaging are suitable to approximate the effective heat conductivity of layered media, as shown in Bringedal and Kumar (2017), they are not applicable to more complex pore structures as also the pore and grain shapes themselves, for example through surface area, affect the effective heat conductivity. We aim at capturing pore-scale effects by incorporating them into a model at the Darcy scale. This is done by deriving the effective heat conductivity. Two of the most common approaches are the method of volume averaging, and the theory of homogenization. When using volume averaging, effective quantities are obtained based on closure relations (Whitaker 1999). Applications to heat transport in porous media can be found in Hsu (1999) and Quintard et al. (1997). However, in the present work we consider the approach of formal homogenization to obtain effective heat conductivities (Auriault et al. 2009;Hornung 1997). Formal homogenization allows to derive upscaled equations and corresponding effective quantities following suitable assumptions. Within this framework, the effective properties of the porous medium are obtained by solving so-called cell problems at the pore scale. If a porous medium exhibits local periodicity in its structure, we can obtain accurate results for general grain shapes using the theory of homogenization (Auriault 1983). Then, representative sections of the porous medium need to be considered at the pore scale and are the basis for deriving upscaled equations with corresponding effective quantities.
In the following, thin porous media are considered. With a thin porous medium, we mean a porous medium where the thickness of the porous medium is of the same order of magnitude as the length of a representative elementary volume at the pore scale. This corresponds to the proportionally thin porous medium (PTPM) considered in Fabricius et al. (2016). Note that different types of thin porous media are also discussed in that paper, in the context of fluid flow. For thin porous media, the boundary conditions applied at the top and bottom boundaries will have an influence on the overall effective behavior and hence on the effective quantities of the porous medium.
The presence of thin layers embedded in a surrounding porous medium and their effect on effective quantities have been analyzed by homogenization in the context of diffusion-reaction systems (Bhattacharya et al. 2020;Gahn et al. 2021Gahn et al. , 2017Gahn et al. , 2018Neuss-Radu and Jäger 2007) and for unsaturated flow List et al. 2020). In these works, the focus has been in particular on transmission conditions across the thin layer and the interaction with the surrounding porous medium as the width of the layer approaches zero. Homogenization techniques can also be applied to free-standing thin structures. By considering a thin strip as a simplified representation of a porous medium, effective equations have been found by combining asymptotic expansions and transversal averaging for reactive transport (van Noorden 2009), heat transport (Bringedal et al. 2015), biofilm growth (van Noorden et al. 2010) and two-phase flow (Lunowa et al. 2021;Sharmin et al. 2020). However, in these works the thin strip is in practice a channel and does not contain a porous structure. On the other hand, when considering a thin porous medium, both the porous structure as well as the boundary conditions at top and bottom boundaries influence the effective behavior.
Much work has been done considering flow in thin porous media. Laminar as well as turbulent single-phase flow in thin porous media have been investigated for example in Chen and Papathanasiou (2008), Fabricius et al. (2016), Hellström et al. (2010), Koch and Ladd (1997) and Wagner et al. (2021). A benchmark comparison for a homogenization-based approach applied to flow through thin porous media in Wagner et al. (2021) shows a good agreement between the results of three-dimensional homogenization and pore-scale resolved models as well as experimental findings. In this paper, we extend the formal homogenization approach applied to fluid flow in thin porous media in Fabricius et al. (2016) by incorporating heat transport and focus in particular on the role of the arising effective heat conductivity.
This work is also related to Bringedal and Kumar (2017), which takes a similar approach to investigate the influence of chemical reactions on effective heat conductivities in twodimensional periodic porous media. In Bringedal and Kumar (2017), the focus is, however, on the impact of the evolving pore structure on the effective quantities. Here, we will not explicitly account for reactive transport and subsequent evolving pore geometries, but we will investigate the effect of different grain shapes and sizes. Moreover, we consider a fully three-dimensional approach. As a consequence, since the porous medium is thin, the type of boundary conditions applied in the third spatial dimension will impact the results. Although we will not explicitly account for any evolution of the grains in this paper, the resulting models for effective heat conductivity found in this paper would under suitable assumptions be the same if allowing the grains to evolve, and complementing the model by an appropriate evolution equation. As shown in Bringedal et al. (2016), the influence of the evolving grains appears in the derived Darcy-scale model, while the effective quantities only depend on the grain shape at a given time (Bringedal et al. 2016;Bringedal and Kumar 2017). Hence, extensions to porous media evolving due to mineral precipitation and dissolution reactions can be formulated and included in the upscaling procedure by considering, e.g., a level-set formulation as in Bringedal et al. (2016).
In the following, we first introduce the model formulation for a thin, periodic porous medium in Sect. 2. The two-scale asymptotic expansions of the formal homogenization approach are introduced in Sect. 3. In Sect. 4, we obtain the respective cell problems at the pore scale. Upscaled equations valid at the Darcy scale are derived in Sect. 5. The Darcyscale equations rely on effective quantities, which are found through solving the corresponding cell problem. We continue with a discussion on the impact of different boundary conditions and grain shapes on the derived effective heat conductivity as well as the impact of the upscaling procedure itself in Sect. 6, before we end with some concluding remarks.

Structure of Porous Medium
The modeled porous medium consists of void space f that is occupied by fluid, and grain space g such that = f ∪ g ∪ S , where S denotes the internal boundaries between void and grain space. The porous medium is characterized by a coordinate system (x 1 , x 2 , z) , and we assume that x 1 , x 2 ∈ [0, L] and z ∈ [0, H] , with H ≪ L since we are dealing with a proportionally thin porous medium. We introduce the parameter to separate Darcy and pore scale. The quantity l describes the horizontal extent of the chosen pore-scale reference domain and defines the length size of the pore scale. Hence, we can model the periodically oscillating domain at the pore scale by introducing the local variable For a given point (x, z) ∈ , we then have (y, z) ∈ Q∶ = [0, l] 2 × [0, H] . We use the superscript to stress that the domain contains highly oscillatory structures. A so-called cell Q, as shown in Fig. 1, consists of void space Q f , which is occupied by fluid, and grain space Q g . The internal surface between Q f and Q g is denoted with S Q . Furthermore, we introduce Q for the internal sidewalls as well as Q for top and bottom of the cell as visualized in Fig. 2. We have that Q = Q f ∪ Q g , where Q f and Q g are the fluid and grain part of the top and bottom surfaces, respectively. Hence, we define the boundary of the void space as Q f = Q f ∪ S Q f ∪ Q and of the grain space as Q g = Q g ∪ S Q g . For simplicity, we here assume that the internal sidewalls Q are occupied by fluid, but allowing grains to touch these sidewalls represents no practical difference for the presented model as long the void space remains connected. Note that S Q f and S Q g indicate the same surface but with opposite orientation. Altogether, we have the surface of the whole cell to be given as Q = Q ∪ Q . We consider Q to be a zoomed-in view for a given x and in Sect. 3 we assume that functions defined on Q are periodic in y . Hence, is built up through several cells Q. Note that these cells do not need to be identical as, e.g., grain shape could vary, but all cells follow the structure as described above. Different grain shapes correspond to the pore-scale geometry varying with x . However, since we will rely on local periodicity in y , two nearby-lying cells Q should not be too different from each other. y = 1 x.

Fig. 1
Darcy and pore scale with cell Q. Note that we consider the z-direction as the vertical direction, but is for simplicity depicted as pointing out of the paper in figure The whole porous medium is defined as the union of all pore-scale cells Q, meaning = ⋃ Q . Accordingly, we obtain f = ⋃ Q f and g = ⋃ Q g . Similarly, the union of all boundaries is defined by

Equation Set and Boundary Conditions
To model fluid flow and heat transport on the pore scale of the porous medium, we formulate corresponding conservation equations and boundary conditions. We consider the conservation of mass and the conservation of momentum by using the Stokes equations where f (T f ) denotes the density of the fluid, u the velocity, p the pressure and (T f ) the viscosity of the fluid. The fluid density and viscosity are assumed to depend smoothly on fluid temperature. We use here the Stokes equations instead of the more general Navier-Stokes equations as we in the upscaling only consider the creeping flow regime to be in the range of Darcy's law. However, with suitable assumptions on the Reynolds number, Darcy's law can also be derived from the Navier-Stokes equations (see, e.g., Bringedal et al. 2016).
To describe the conservation of energy in the system, we consider a simplified system. Since we are considering creeping flow in a porous medium, fluid velocities will generally be low. Further, we consider temperature ranges such that no phase change occurs. In this case, energy conservation can be formulated in terms of temperature (Landau and Lifshitz 1987). We introduce the two temperatures T f defined in the void space f , and T g defined in the grain space g . In the fluid, heat transport is due to advection and conduction, while any internal heat production due to friction is neglected. The heat transport in the grain space only occurs through conduction: The quantities k f and k g denote the heat conductivities of the fluid and the grain, whereas c f and c g are the specific heats. They are all assumed to be material constants. Finally, g (T g ) is the density of the grain, which is assumed to depend smoothly on the grain temperature. Note that the equations (1-4) are strongly coupled as the fluid flow and temperatures jointly affect each other through the advective term, and since the densities and viscosity are temperature dependent.
Note that although we here consider relatively simple energy equations (3), (4), extending the model by including, e.g., temperature-dependent specific heats in the upscaling is straightforward as long as they depend smoothly on temperature. However, due to the special role of heat conduction, allowing the heat conductivities to depend on temperature would affect the upscaling procedure as higher-order effects for this process play a role. Hence, we are limited to consider moderate temperature ranges and fluids and solids such that the assumptions behind (3), (4) remain valid.
To close the model, boundary conditions at all boundaries as well as initial conditions for the variables are needed. Since the upscaling only depends on boundary conditions on and S , only these will be specified in the following. Regarding the flow, we assume no-slip boundary conditions on the grain surface and on the top and bottom boundaries of the porous medium: No-slip conditions are commonly used for porous-media flows and ensure that the fluid flow remains in the creeping regime. For completeness, we mention that extensions with slip also exist (Lasseux and Valdés-Parada 2017;Lunowa et al. 2021). The impact of Neumann as well as Dirichlet boundary conditions on the top and bottom boundaries of the domain for the energy conservation will be considered: Hence, we will use either (6a) and (7a), or (6b) and (7b). Here, e z denotes the basis vector in z-direction. Further, we have introduced the heat flux and the temperature that are specified on top and bottom boundaries of the cell Q. The Neumann boundary conditions correspond to an applied heat flux at the top and bottom boundaries. Applying Dirichlet boundary conditions corresponds to assuming that the top and bottom boundaries are perfectly heat conducting and hence adopts the applied temperature imposed here. These two types of boundary conditions represent the two complementary cases of the top and bottom boundaries either transmitting a given heat flux or being perfectly heat conducting. Combinations of these two cases can be formulated through Robin boundary conditions, but are outside the scope of this work. Note that the heat flux or the temperature could themselves come from a heat transport model for an adjacent medium, which is why they could in the general case vary with time and space. We are interested in how the two different types of coupling to the thin porous medium affect the heat transport. Note that the boundary conditions are for consistency assumed to be the same on the top and bottom boundaries. Furthermore, they are assumed to not depend on the local variable y , meaning that they cannot be highly oscillating in space. We also assume that they are not highly oscillating with time. The latter two assumptions are needed to enable a separation of scales when upscaling the model.
On the interface between fluid and grain, we assume local thermodynamic equilibrium: which is a reasonable assumption in the creeping flow regime. Alternatively one can apply a contact conductivity model, and we refer to Auriault et al. (2009)[Chapter 4.3] for the influence of such boundary conditions. To conserve the energy, the heat flux from fluid to grain is the same as from grain to fluid: where n denotes the normal vector on S .

Non-Dimensionalization
The equations and boundary conditions (1) -(9) are non-dimensionalized by introducing the following reference quantities to define non-dimensional variables and quantities: Since we consider thin porous media, we assume that H, which is the length of the domain in z-direction, is of the same order as the length l of the previously defined cell. This leads to H ∶= H ∕l = O( 0 ) , where the non-dimensional height H does not depend on . We use the horizontal length scale l to non-dimensionalize the vertical direction to emphasize their same order in size. Our assumption corresponds to the case 'Proportionally Thin Porous Medium' introduced in Fabricius et al. (2016). The non-dimensionalization of the spatial variables in (10) allows us to introduce the non-dimensional domains Due to the different scalings used for x and ỹ , both pore-scale and Darcy-scale domains appear as unit sized domains. Note that ∇ denotes the non-dimensionalized nabla operator in the following. By inserting the non-dimensionalized definitions (10), all model equations and boundary conditions can be rewritten in terms of non-dimensional variables and characteristic non-dimensional numbers. From dimensional analysis, we expect 5 independent nondimensional numbers in addition to the length scale ratio . These 5 non-dimensional numbers are described and chosen below. The conservation of mass (1) and momentum (2 Since we use Stokes equations, it is implicitly assumed that the Reynolds number is small. In order to ensure that we are in the creeping flow regime where Darcy's law is valid, we assume ReEu = O( −2 ) . Hence, we write ReEu = k flow −2 , where k flow is a non-dimensional constant not depending on . The equations describing the heat transport in the void (3) as well as the grain space (4) are rewritten in the following way (11) t̃f +∇ ⋅ (ũ̃f) = 0 iñf , Here, Pe f and Pe g are the thermal Péclet numbers for fluid and grain, and hence, our last two independent non-dimensional numbers. We assume Pe f , Pe g = O( 0 ) , which means that the time scales for heat conduction and heat advection are the same. Our assumption implies that heat transport by conduction and advection are equally important at the macroscopic scale, and both will appear in the resulting upscaled model. Different assumptions on the Péclet number would lead to either only conduction (low Péclet number) or only advection (high Péclet number) to be dominating at the macro-scale. Note that our assumption also implies that the ratio between k f and c f has to be the same order of magnitude as the ratio between k g and c g . For the influence of large differences between the heat conductivities, we refer to Auriault et al. (2009) [Chapter 4.2]. In the following, we denote with f and g the non-dimensional heat conductivities of fluid and grain, respectively.
The non-dimensionalized boundary conditions are: Since we will only consider the non-dimensional variables in the following, we skip the tilde.

Two-Scale Asymptotic Expansions
To incorporate the dependence on the local variable y , we introduce two-scale asymptotic expansions for the velocity, the pressure as well as the temperatures in the fluid and the grain.
We hence have that these functions i are periodic in y , meaning that for all (y 1 , y 2 , z) ∈ Q . Due to the difference in scaling for x , y and z, the non-dimensional nabla operator is where For the non-dimensional Laplace operator Furthermore, we introduce for convenience the short-hand notation Recall that the densities and the viscosity are assumed to depend smoothly on the temperature. By inserting the asymptotic expansions for the temperatures and applying the Taylor expansion, we get where we have introduced the short-hand notations f0 , f1 , g0 , g1 , 0 and 1 to describe the temperature-dependence as indicated in (20)-(22). We are interested in the limit as → 0 and analyze which terms of the model equations are then dominating. By inserting the asymptotic expansions into the model equations and sorting the terms after increasing order in terms of , we can isolate those terms.
When the asymptotic expansions of the velocity and fluid temperature are inserted into the conservation of mass (11), the terms for the two lowest orders in are: (20) Similarly, the conservation of momentum (12) results in For the conservation of energy within the void space (13), we insert asymptotic expansions for the velocity as well as the fluid temperature and get Using the respective asymptotic expansions in the non-dimensionalized conservation of energy in the grain space (14), yields The two lowest orders of the velocity's no-slip boundary are obtained by inserting the asymptotic expansions of the velocity into the non-dimensionalized boundary condition (15): Analogously, the Neumann and Dirichlet boundary conditions for the temperature on the top and bottom boundary of the void and grain space (18) and (19) result in

3
where we for convenience have multiplied the Dirichlet conditions (18b) and (19b) with −1 . The internal boundary condition for the temperature in void and grain space (16) holds in a similar way for the different orders as well: By inserting the asymptotic expansions for the temperatures into the internal flux boundary (17), we get

Dominating Order Behavior
Based on the dominating term in the momentum equation (25), we can conclude that which means that p 0 does not depend on y and z. Equivalent results for T f0 and T g0 follow from the dominating terms in the conservation of energy (27) and (30), and the respective top and bottom boundaries (35a) -(36b) as well as the interface condition (43). Hence, T f0 and T g0 are constant within a cell and fulfill because of the continuity condition on the interface (41). Moreover, since 0 , f0 and g0 only depend on T f0 and T g0 , respectively, this leads to Hence, these variables remain constant within the respective cell for a given x at the Darcy scale.

Cell Problems
The equations corresponding to a fixed order with respect to from the previous section are now used to derive so-called cell problems which describe the variations at the pore scale. We first consider the conservation of mass and momentum which results in a cell problem for the flow. The corresponding derivation for the upscaled flow in a porous medium can be found, e.g., in section 1.4 of Hornung (1997) or for the proportionally thin porous medium in Fabricius et al. (2016), and is only included for completeness here. Our main focus is on the conservation of energy and the corresponding cell problems for heat transport which are discussed afterwards.

Conservation of Mass and Momentum
Due to the linearity of the terms of order 0 in the conservation of momentum (26), we can rewrite p 1 = p 1 (x, y, z, t) and u 0 = u 0 (x, y, z, t) as linear combinations of x j p 0 for some unknown weights r j and s j . In order to ensure the periodicity of p 1 and u 0 , we require periodicity of r j and s j in y.
We insert the linear combinations (49) and (50) into the second lowest order terms of the conservation of momentum (26) as well as the lowest order terms of mass conservation (23) and the second lowest order terms of the no-slip boundaries (34). By making use of the dependencies of density and viscosity in (48), we obtain the following cell problems to determine r j and s j for j = 1, 2: where e j denotes the unit vector in x j -direction. As mentioned before, s j and r j are periodic on Q . We add the constraint to ensure the uniqueness of the solution. One should note that this cell problem only needs to be solved in y and z. The effect of x comes in as a parameter accounting for the chosen pore geometry at a position x at the Darcy scale. Hence, the cell problem accounts for the local flow at the pore scale and will later be used to define the permeability of the porous medium. By solving for r j and s j in ( P flow (x) ), we can obtain p 1 and u 0 as functions of ∇ x p 0 through (49) and (50).

Conservation of Energy
Next, we apply a similar approach to the conservation of energy: The linearity of the terms (28) and (31) of order −1 indicates that T f1 and T g1 can be represented as linear combinations of x j T 0 , factorizing the dependencies on x and y . Since the following steps depend on which boundary conditions are used on the top and bottom boundary Q , we divide our approach into two cases: We first consider Neumann boundary conditions, and then, the results for Dirichlet boundary conditions are presented. For a clear distinction between the two cases, we use the superscript 'N' or 'D', respectively, to indicate the considered boundary conditions.

Neumann Boundary Conditions on Top and Bottom Boundary Q
Due to the linearity of the problem, we choose the following expressions in order to satisfy the boundary conditions at the top and bottom boundaries: where v N j and w N j ( j = 1, 2 ) are weights that are to be determined. Note that the last term in (51a) and (51b) ensures that the Neumann boundary conditions (37a) and (38a) are fulfilled. In the case of homogeneous Neumann boundary conditions, the last term is zero. By inserting the linear combination (51a) into the second lowest order terms of the energy conservation (28) and using the dependencies of T f0 in (47) as well as the lowest order terms of the mass conservation (23), we get For the grain space, a similar expression can be derived when the linear combination (51b) and the respective equation for energy conservation (31) are used: Suitable boundary conditions result from inserting the linear combinations into the second lowest order terms of the existing boundary conditions (37a), (38a), (41) and (44): By combining the derived equations (52) -(57), we get two coupled cell problems, for j = 1, 2: The weights v N j and w N j are required to be periodic in y due to the periodicity of T f1 and T g1 . We add the uniqueness constraint with an arbitrary constant. Note that since we added the inhomogeneous Neumann boundary condition as part of the series expansions (51), we obtain a cell problem which is independent of the heat flux applied at the top and bottom boundaries. Hence, ( P heat, N (x) ) is used both for homogeneous and inhomogeneous Neumann boundary conditions. Dimensionality of the Cell Problem If we consider a pore geometry that is constant with respect to z, we find that which means that we can reduce the three-dimensional cell problem to two spatial dimensions and still obtain the same results. To see this, we start by rewriting where we have omitted writing the dependence on x to simplify the notation. Consequently, we can rewrite the cell problem, namely the Laplace equation in Q f (52) as well as the Neumann boundary condition (54), based on the separation of the variables y and z, as Hence, we have derived two equations that involve either only y or z: There exist only two types of solutions for o j f (z) that fulfill the z-dependent equation (63) as well as the Neumann boundary condition in (62): Note that since (65) leads to a constant solution in case of c = 0 , we only consider c ≠ 0 here. Using same arguments for Q g and Q g , we find that o j g (z) has a form corresponding to (64) or (65). In case of a constant solution (64), it is trivial that this solution does not depend on z and satisfies (62). Hence, we now show that (65) cannot be a solution of our problem.
As mentioned above, v N j and w N j need to satisfy a uniqueness condition as given in (58) for an arbitrary constant. In the following, we set the constant to one. If we incorporate the separation approach (61), we get However, the first term can be rewritten into where Q f,2D is a two-dimensional cross section of the fluid domain in the cell. Note that this separation can only be done since the geometry does not vary with z. Similarly, we rewrite the integral of w N j . The uniqueness condition is calculated as which is a contradiction. Hence, o j f (z) as well as o j g (z) have to be constant functions as stated in (64). The three-dimensional cell problem can therefore be reduced to a two-dimensional problem if the geometry does not depend on z. Note that we still obtain three-dimensional solutions for the cell problem that cannot be reduced to two-dimensional cell problems, if the grain shape changes along the z-axis.

Dirichlet Boundary Conditions on Top and Bottom Boundary Q
To satisfy the Dirichlet boundary conditions on Q , we can use a similar linear combination as for zero-Neumann boundary conditions: where v D j and w D j ( j = 1, 2 ) are weights that are to be determined. We insert the linear combinations in the respective equations following similar steps as before and obtain the following coupled cell problem for j = 1, 2: This cell problem differs from the previous case of Neumann boundary conditions in the boundary conditions on Q f and Q g . Note that due to the Dirichlet boundary conditions, these cell problems cannot be reduced to two-dimensional even if the geometry is not depending on z.

Darcy-Scale Equations
In the following, we integrate the flow equation (50) over its domain Q f . Further, the next order terms with respect to from the mass (24) and energy conservation (29) and (32) are integrated over the respective domains within Q. This is done in order to obtain equations at the Darcy scale, while including effects from the pore scale through effective quantities.

Effective Flow
The derivation of the effective flow, meaning Darcy's law, can be found in the literature (e.g., in Sect. 1.4 of Hornung (1997) in the case of constant viscosity and density, and Bringedal et al. (2016) with temperature-dependent viscosity and density) and is again only shown for completeness. For the effective flow ū 0 in the porous medium, the definition of u 0 as stated in (50) is integrated where the components of K are and the viscosity follows the dependency in (48). The components s j,i follow from solving the local, three-dimensional cell problem ( P flow (x) ). The averaged equation (68) is Darcy's law, where K is the permeability of the porous medium. The matrix K is symmetric and positive definite and we refer to Lemma 4.2 in Sect. 1.4 in Hornung (1997) for a detailed proof. Note that the permeability in our case is a 2 × 2 matrix, although the cell problem is solved in a three-dimensional domain. In Fabricius et al. (2016) and Wagner et al. (2021), the authors discuss the applicability of approximate, two-dimensional forms of the cell problem, accounting for the porous medium being thin. Hence, under suitable assumptions, the resulting permeability for thin porous media can be found more cheaply.
By integrating the second lowest order terms of the mass conservation (24) over the volume and using Gauss' theorem, the following equation is derived: Note that the boundary of the void space was previously defined as The integral over Q is zero due to the periodicity of u 0 and u 1 . In addition to that, the integral over Q f and S Q f disappears because we assumed no-slip boundaries in (34). The upscaled equation (70) can be rewritten based on the definition of ū 0 in (68) as where ∶= |Q f | ∕|Q| is the porosity.

Effective Energy Conservation
To describe the conservation of energy in the void as well as the grain space at the Darcy scale, the third lowest order terms in void and grain space (29) and (32) are used. If we integrate them over their respective domains, add them up, divide by |Q| and apply Gauss' theorem as well as the definition of the effective flow (68), we obtain Considering the composition of Q f = Q f ∪ S Q f ∪ Q and Q g = Q g ∪ S Q g , we first show that term C vanishes: For the integrals over Q f ∪ S Q f , we make use of the no-slip boundary conditions (33) and (34) again. Since all functions within the integral are periodic in y , the contributions of the different parts of Q cancel each other.
Note that Q f ∩ Q g = S Q , while the respective normal vectors on this domain point in opposite directions for the fluid ( S Q f ) and the grain component ( S Q g ). Due to the continuity condition (45), the contributions from S Q in term B cancel each other. The integral over Q is zero due to periodicity. Hence, term B results in where we have introduced n 3 = ±1 to denote the third component of the normal vector on Q f and Q g , respectively. The following steps depend on the top and bottom boundary of the domain. Hence, we consider the case of Neumann boundary conditions first before moving on to Dirichlet boundary conditions.

Neumann Boundary Conditions on Top and Bottom Boundary Q
In this case, the contributions of the integrals over Q f and Q g in term B as given in (73) disappear according to the Neumann boundary conditions (39a) and (40a). Therefore, term B is zero.
To simplify the terms A 1 and A 2 , we insert the linear combinations (51a) for T f1 and (51b) for T g1 and use that T f0 as well as T g0 do not depend on y and z:

with components
We can further simplify these expressions by using the following relation: This relation is obtained by rewriting the left-hand side of the equation and applying the Gauss' theorem: The integral over Q vanishes due to periodicity of v N j , and on Q f the normal component n i is zero since the respective normal vector is parallel to the z-axis. If we use the continuity condition (56) of the cell problem, we get

3
because S Q f and S Q g denote the same surface but with opposite orientation. When applying the previous steps in inverse order, we end up with Hence, we have (76). Making use of the relation between the integrals (76) and the quantities V N and W N in (74) and (75), the upscaled effective heat conductivity is given as with components i, j = 1, 2 Note that the first term represents the volume-weighted average of the heat conductivities in fluid and grain. The second term accounts for the internal structure of the heat transport within the cell problem and the interaction with the top and bottom boundaries, which are not accounted for by the volume-weighted averages.
If we insert all derived expressions above into (72), we obtain the upscaled form of conservation of energy:

Dirichlet Boundary Conditions on Top and Bottom Boundary Q
We first consider the simplified form of term B as stated in (73): If we assume that T f2 and T g2 are symmetric with respect to z, term B is zero since the remaining integrals sum up to zero. This assumption is reasonable due to the symmetric boundary conditions and since the setup of the discussed problem yields a symmetric behavior.
As done for the case of Neumann boundary conditions, we also insert the respective linear combinations (67a) and (67b) into the terms A 1 and A 2 . This results in: The components of V D and W D are given by: Note that the integrals of y i v D j over Q f and of y i w D j over Q g satisfy The relation can be derived analogously to the case of Neumann boundary conditions (see (76)). By using this relation as well as the definitions of V D and W D , we define the effective heat conductivity with components i, j = 1, 2 As in (77), we see that the effective heat conductivity can be written as the sum of the volume-weighted average of the heat conductivities, and of a part accounting for the internal structure of the heat transport and interactions with top and bottom boundaries. Hence, we obtain the following Darcy-scale equation for the heat transport:

Remarks Regarding the Effective Heat Conductivities
If we compare the effective heat conductivities S N in (77) for Neumann boundary conditions and S D in (79) for Dirichlet boundary conditions on the top and bottom boundary, one finds that they have the same structure. However, they depend on the solution of the respective cell problems ( P heat, N (x) ) and ( P heat, D (x) ), which are not the same. Hence, we obtain different effective heat conductivities depending on the type of boundary conditions. The effective heat conductivities S N in (77) and S D in (79) are both symmetric and positive definite, as one also finds for two-or three-dimensional porous media where the top and bottom boundary conditions do not have an influence Auriault (1983). To see this, one considers the weak form of the cell problems ( P heat, N (x) ) and ( P heat, D (x) ). Using test functions that also fulfill the top and bottom boundary conditions, the weak form for ( P heat, N (x) ) is for all sufficiently smooth that are periodic in y and fulfilling z = 0 on Q . Choosing equal to v N j in Q f and to w N j in Q g , and using (76) together with Gauss' theorem, we obtain This identity can be used to rewrite the components S N ij to .
which is obviously symmetric, and can be seen to be positive definite by considering the sum ∑ 2 i,j=1 i S N ij j for real numbers i (see also Proposition 3.2 for a diffusion problem in Hornung 1997). For the cell problem with Dirichlet boundary conditions ( P heat, D (x) ), the argument follows similar steps, but using test functions fulfilling = 0 on Q . The corresponding weak form and components of S D ij can be written the same way as for the Neumann boundary conditions case, Note that although the cell problems ( P heat, N (x) ) and ( P heat, D (x) ) generally need to be solved in three-dimensional domains, the resulting effective matrices S N (x) and S D (x) only vary in the horizontal directions, meaning along x = (x 1 , x 2 ) . A z-dependence of the effective heat conductivity would appear if higher order terms from the conservation of energy in (13) and (14) would be included, which would correspond to a better approximation of the original problem. The presented effective model is based on terms up to the order O( 0 ) . In the expressions for A 1 + A 2 , the divergence with respect to x is zero in the third component. Therefore, the third row and column of the effective heat conductivity are zero and we only need to consider S ∈ ℝ 2×2 .
Hence, in a thin porous medium, the resulting effective heat conductivity for the horizontal heat transfer depends on the three-dimensional structure of the pore scale and on the boundary conditions applied to the top and bottom boundaries of the thin porous medium.

Summary of Upscaled Model
The derived upscaled model consists of Darcy's law (68) and upscaled mass conservation (71), together with an upscaled equation for the effective heat transport (78)  where S = S N in the case of Neumann boundary conditions on the top and bottom boundaries, and S = S D in the case of Dirichlet boundary conditions. Equations of state describing how the fluid and grain densities and fluid viscosity depend on temperature must also be included to close the system. Note that these equations are all defined on two-dimensional domains, as the vertical (thin) direction does not need to be resolved. Hence, the effective matrices K, S are 2 × 2 matrices. The components of the permeability matrix K are determined by (69), where the solution of the corresponding cell problems is given through ( P flow (x) ). The components of the effective heat conductivity S N and S D are given by (77) and (79), which depend on the corresponding cell problems ( P heat, N (x) ) and ( P heat, D (x) ), respectively. Note that the cell problems are defined on three-dimensional domains, which are small portions of the original domain. Despite the strong coupling between the original model equations, the cell problems can be solved independently of each other.
The upscaled model is valid under the assumptions on the original pore-scale model equations stated in Sect. 2.2, and under the assumption of the non-dimensional numbers as explained in Sect. 2.3. Note that the (scaled) non-dimensional numbers still appear in the resulting upscaled model through the effective quantities K, S and through the cell problems ( P heat, N (x) ) and ( P heat, D (x)).

Effective Heat Conductivity Behavior Based on Cell Problems
In the following, we analyze the effective heat conductivities that are calculated based on the solutions of the cell problems ( P heat, N (x) ) and ( P heat, D (x) ). For the numerical results below, we have used Netgen for mesh generation and the finite element software NGSolve to solve the weak form of the cell problems (see Schöberl 2014). Note that the cell problems ( P heat, N (x) ) and ( P heat, D (x) ) are elliptic. To generate the mesh in Netgen, we specify the domain and location of the inner boundary and give the maximum mesh size h max . However, Netgen generally employs smaller grid cells near inner boundaries. To discretize the cell problems on these meshes, we use subspaces of H 1 , using polynomials up to third order as basis functions. The remaining computations to obtain the effective heat conductivities according to (77) and (79) have been performed using NGSolve, as the needed solution derivatives are directly available through NGSolve. As mentioned above, the cell problem can be reduced to two dimensions whenever we apply Neumann boundary conditions on the top and bottom boundary Q and consider a setup where the grain shape does not change along the z-axis. By reducing h max , larger accuracy is obtained. We found that four significant digits in the cell problem solution and in the effective quantity was obtained already by h max = 0.1 . The number of used grid cells was in this case for the cell problems in Sect. 6.1-6.3 in the range 180-400 for two-dimensional cell problems, and 6 000-144 000 for three-dimensional cell problems. The largest amount of grid cell were needed for ellipsoid-shaped grains.
The provided code Scholz and Bringedal (2021) calculates all components of the effective heat conductivities S N and S D . We will limit our attention to isotropic grain shapes. In this case, the off-diagonal components S N ij and S D ij ( i ≠ j ) are close to zero and S N 11 = S N 22 as well as S D 11 = S D 22 . Therefore, we only present and discuss the first diagonal component in the following.
We are interested in which way the use of Neumann or Dirichlet boundary conditions on the top and bottom boundary Q as well as different grain shapes affect the effective heat conductivity. Besides that, different ratios of the heat conductivity of the grain g and the fluid f count among the parameters of interest. To simplify the comparison, we set f = 1 and consider different values for k = g∕ f = g . Note that if follows from assumptions on Pe f and Pe g that f and g are the same order of magnitude. Hence, with f = 1 , g have to be close to 1 as well. In the following, we consider k between 0.2 and 5.

Effect of Boundary Conditions on Top and Bottom Boundary Q
As shown in Sections 4.3 -4.4, the cell problems for effective heat conductivities in case of Neumann and Dirichlet boundary conditions on Q differ. Despite the similar structure of the equations as given in (77) and (79), the respective cell problems ( P heat, N (x) ) and ( P heat, D (x) ) apply different boundary conditions. In the following, we consider circular cylinder-shaped grains of different radii and compare the results for the two boundary condition types on Q . The resulting effective heat conductivities are shown in Fig. 3a.
The results for Neumann boundary conditions are consistent with the ones presented and discussed in Bringedal and Kumar (2017), where two-dimensional cell problems without the effect of boundary conditions were considered. However, this is due to fact that the grain shapes do not vary along the vertical axis in this case. For Dirichlet boundary conditions, we obtain a similar behavior with respect to varying the size of the grain, but overall larger effective heat conductivities. In the trivial case of k = 1 , fluid and grain have the same heat conductivity and we obtain the same constant effective heat conductivity for Neumann and Dirichlet boundary conditions. However, for k ≠ 1 the effective heat conductivity is always larger in case of Dirichlet boundary conditions than for Neumann boundary conditions on Q . This positive impact of Dirichlet boundary conditions is due to them representing perfect heat conduction on Q . If we consider the respective difference in detail for a cylinder and a cuboid, as shown in Fig. 3b, we observe that for → 0 and → 1 , meaning the grain fills almost all of the void space or its size is negligible, this difference goes down to zero. This is the case because the local domain is almost homogeneous and therefore equally heat conductive, independently of the boundary conditions. For porosities that are neither close to zero nor close to one the choice of Dirichlet boundary conditions has a strong, positive effect on the effective heat conductivity.
Note that the surface area between fluid and grain is proportional to the radius when considering cylinder-shaped grain. Hence, Fig. 3a can be read as effective heat conductivity as a function of fluid-grain surface area when scaling the horizontal axis with a factor 2 H = 2 for this setup. Hence, when one knows that the grains are cylinder-shaped, the effective heat conductivity can be determined by the cylinder radius, correspondingly the fluid-grain surface area.
Since effective heat conductivities for porous media are often calculated in terms of porosity-weighted averages of the individual heat conductivities found in fluid and grain (see, e.g., Nield and Bejan 2017), we compare in Fig. 3c the effective heat conductivities found through the cell problems with such averages. The averages are calculated through As seen in Fig. 3c, and as pointed out earlier in Bringedal and Kumar (2017) for twodimensional porous media without the influence of top and bottom boundary conditions, these porosity-weighted averages offer approximate values for the effective heat conductivity, but cannot predict the detailed behavior as the porosity varies.

Effect of Different Grain Shapes
We are interested in understanding to which extent the effective heat conductivity is affected by changes within the cross-sectional shape of the grain along the z-axis. The resulting equations for the effective heat conductivities (77) and (79) directly indicate the impact of the porosity on the effective heat conductivity. Therefore, we compare different grain shapes but always using the same grain volume |Q ref g | in order to account for effects of the detailed shape. Since we consider the same porosity, any porosity-weighted average would always give the same value (for fixed f and g ) independent of the shape. All considered shapes are rotationally symmetric with respect to the z-axis. Hence, we can introduce r(z) ∶ [0, 1] → [0, 0.5] to describe the grain radius perpendicular to the axis for different values of z. Further, the minimum and maximum radius of a grain are given by As grain reference volume |Q ref g | , we use a cylinder with radius r cyl (z) = 0.25 . In addition to the cylinder, we construct a cone-shaped as well as an ellipsoid-shaped grain which have the same volume as our reference grain. They are defined by r min = min z∈ [0,1] r(z), r max = max z∈ [0,1] r(z).
In order to ensure that all grains are of the same volume |Q ref g | , we obtain the following restrictions on the minimum and maximum radii of the cone and the ellipsoid for r min∕r max → 0: In case of the limit r min∕r max → 1 , the shapes of ellipsoid and cone approach the shape of the reference cylinder. The detailed setup is sketched in Fig. 4. The minimum and maximum radius of the different shapes satisfy r con min ≤ r ell min ≤ r cyl and r cyl ≤ r ell max ≤ r con max for any given ratio of the radii.
In the following, we compare the effective heat conductivities for different conductivity ratios k and Neumann as well as Dirichlet boundary conditions on the top and bottom boundary Q . In Fig. 5, we show how the results vary with the ratio r min∕r max , the fluid-grain interfacial area and the contact area between top and bottom boundary and the more conductive phase. We do not show how the results vary with porosity as in Fig. 3, since the porosity is kept constant. Note that since only one cylinder is considered here, the results for the cylinder are included either as a reference line or as a reference point.
Neumann Boundary Conditions on Top and Bottom Boundary Q The dependence of the effective heat conductivity on the grain shape and the corresponding maximum and minimum radius is clearly visible: A variation of the grain shape causes larger effective heat conductivities. For k < 1 (Fig. 5a), the grain is less conductive than the fluid and therefore hinders the overall conduction of heat. Hence, decreasing the minimum radius of the grain causes a better connectivity of the fluid which increases the effective heat conductivity. Recall that in case of an ellipsoid-or cone-shaped grain, the minimum radius is smaller than for the reference cylinder. That is why we observe significant changes in the effective heat conductivity for the ellipsoid and the cone compared to the effective heat conductivity of the reference cylinder for small ratios, especially in the case of a cone-shaped grain, where the smaller minimum radius is found. r con (z) = (r con min − r con max )z + r con max , For k > 1 (Fig. 5b), the grain is more conductive than the fluid. Consequently, having better connected grains, meaning a smaller distance between neighboring grains, increases the effective heat conductivity. This is obtained when the maximum radius of the grains Dirichlet (light-colored) boundary conditions on the top and bottom boundary Q . Results for the respective cylinder are given as reference is increased. As mentioned above, the cone and the ellipsoid are characterized by a larger maximum radius compared to the cylinder, which indicates a better connection between the grains in neighboring cells. Therefore, we obtain larger effective heat conductivities also in this case when varying the grain shape, in particular for the cone where a larger maximum radius is found.
The increasing eccentricity of the cone-and ellipsoid-shaped grains corresponds to larger interfacial area between fluid and grain. Hence, for the Neumann boundary conditions, the effective heat conductivities generally increase with increasing surface area between fluid and grain, as seen in Fig. 5c and 5d. However, the interfacial area and value of k alone are not sufficient to determine the effective heat conductivity, and knowledge about the shape (cone or ellipsoid in this case) is also necessary. As seen from Fig. 5e and 5f, although the effective heat conductivity varies with the contact area between top and bottom and the more conductive phase, no specific trend can be determined when using Neumann boundary conditions. Dirichlet Boundary Conditions on Top and Bottom Boundary Q By changing the boundary conditions on Q , we expect the effective heat conductivities to increase due to the positive impact of Dirichlet boundary conditions as discussed in the previous Sect. 6.1. However, the magnitude and overall importance of the positive impact of Dirichlet boundary conditions depends on the surface area of the fluid (for k < 1 ) or the grain (for k > 1 ), respectively, toward Q . In case of k < 1 , the grain hinders the conduction of heat. Hence, the impact of Dirichlet boundary conditions is expected to be larger if the fluid has a larger surface area toward Q . This corresponds to r(z) being smaller on Q . For k > 1 , the opposite holds: A larger surface area of the grains toward Q increases the effective conductivities under Dirichlet boundary conditions. That is achieved by a larger radius r(z).
We observe the same behavior for k < 1 (Fig. 5a) as in the corresponding case of Neumann boundary conditions but with the mentioned general increase based on the different boundary conditions. However, the interplay between increased fluid-grain interface area and the contact area between fluid and top and bottom is not straightforward. In Fig. 5c and 5e, we see how the ellipsoid generally shows an increasing trend with both interfacial area and contact area to top and bottom, while the cone is dominated by the increased interfacial area.
For k > 1 (Fig. 5b), the effective heat conductivities are still always larger than in the corresponding case of Neumann boundary conditions on Q . However, the values for an ellipsoid-shaped grain are now significantly smaller than for the reference cylinder. An ellipsoid-shaped grain has less surface area toward Q , especially for small ratios and therefore small r ell min . In case of a cylinder, the grains exhibit a larger surface area toward Q since r cyl ≥ r ell min . The results for a cone-shaped grain remain above the ones for the cylinder since the cone-shaped grain has a large surface area to the bottom boundary. This is visible in Fig. 5d and 5f. The effective heat conductivities for the cone-and ellipsoid-shaped grains are dominated by the increased contact area between grain and top and bottom boundaries. Here, the ellipsoid's effective heat conductivity therefore decreases with increased fluidgrain interfacial area.
Summary and Implications for Effective Heat Conductivities Both in the case of Dirichlet and Neumann boundary conditions at the top and bottom boundaries, information about the fluid-grain surface area and/or contact area to the top and bottom boundaries and porosity are not sufficient to determine the effective heat conductivitiy. Knowledge of the grain shape is also needed. Depending on the shape, different trends or dependencies on these surface areas are expected. As observed in Sect. 6.1 for circular cylinders, common porosity-weighted averages of the fluid and grain give at best an approximation of the overall effective heat conductivity and cannot account for different shapes nor boundary conditions. Hence, if accurate knowledge of the effective heat conductivity is required, solutions of the cell problems ( P heat, N (x) ) and ( P heat, D (x) ) are necessary.

Comparison Between the Original and Upscaled Model
We here address the influence of the upscaling procedure by comparing the results from the original pore-scale model on an oscillating three-dimensional domain with the upscaled model on an effective two-dimensional domain using the estimated effective heat conductivity. Since we want to focus on the effective heat conductivity, we design a simplified setup employing heat conduction only. Comparison between pore-scale and upscaled approaches for flow can be found, e.g., in Wagner et al. (2021). The domains for the original pore-scale and upscaled models are shown in Fig. 6. The full pore-scale domain consists of 14 × 10 cylinders, with radius 0.2 and 0.4 in the left and right half of the domain, respectively. This domain corresponds to = 0.1.
We consider steady-state heat conduction, hence we solve on the three-dimensional pore-scale domain (Fig. 6a), and on the two-dimensional upscaled domain (Fig. 6b). The temperatures T f and T g from the pore-scale model will give the full variability, including any local oscillations, while the effective temperature T from the upscaled model can only account for the average behavior.
For the left and right boundaries, we apply Dirichlet boundary conditions, using T 1 = 0 and The planes A and B indicate where solutions will be plotted along later (b) Two-dimensional upscaled domain Fig. 6 Computational domains for the pore-scale (left) and upscaled (right) models. The pore-scale domain consists of 14 × 10 cylinders, with two different radii. The upscaled domain hence has two different effective heat conductivities, denoted S (1) 11 and S (2) 11 . The left and right boundaries have Dirichlet boundary conditions, using T 1 = 0 and T 2 = 1 . The other boundaries have homogeneous Neumann boundary conditions T 2 = 1 (see Fig. 6), while the long-sides boundaries have homogeneous Neumann boundary conditions. For the pore-scale domain, the top and bottom boundaries fulfill homogeneous Neumann boundary conditions. We set f = 1 and g = 5 , which is also used to find the effective heat conductivity in the two halves of the upscaled domain through ( P heat, N (x) ) and (77). Since the upscaled heat conductivites are isotropic tensors, we only need the first component and we have that S (1) 11 = 1.183 and S (2) 11 = 2.021 for the left and right halves of the upscaled domain, respectively. Both domains are meshed using Netgen and solved with finite elements through third order subspaces of H 1 using NGSolve. Since h max = 0.1 was found to be sufficient for the accuracy of the cell problems, we employ that here for the full pore-scale domain. This leads to a total of 782 201 grid cells as many internal boundaries need to be resolved by the mesh. The upscaled domain is two-dimensional and meshed with 318 grid cells. By comparing with (upscaled) solutions on finer meshes, we estimate the accuracy of the upscaled model to be five significant digits on this coarse mesh.
Since the temperatures in the pore-scale domain do not vary with the vertical direction (but with both horizontal directions), we can plot the solution along the planes A and B (marked in Fig. 6a) as lines. The effective temperature in the upscaled model does not vary in the x 2 -direction. The full solutions along the planes A and B and the effective solution from the upscaled model along the x 1 -direction are shown in Fig. 7. The temperatures are found to match very well, although small-scale oscillations are found in the full solution which are not captured by the upscaled model. This is expected since we use a nonzero value of . From the difference plot in Fig. 7b, the changes are found to be larger when considering a plane crossing through both fluid and grains in the pore-scale domain (plane B). However, the average of the deviations for each unit cell is close to zero, showing that the upscaled model captures the average behavior almost perfectly. Increased deviations near the external boundaries and near the transition at x 1 = 7 are not observed. Comparing the heat flux in the x 1 -direction (not shown) gives a similar deviation between the full pore-scale solution and the effective solution as shown in Fig. 7b, with maximum deviation around 0.08.
Note in particular that for the same estimated accuracy, the full pore-scale domain needed over 780 000 grid cells, while the upscaled domain was discretized using only 318 grid cells. To find the effective heat conductivities, we have to solve local cell problems as well. For this setup, the cell problems were discretized using 218 and 184 grid cells. Here, (a) Temperature across respective domains from the upscaled model and along two planes in the full pore-scale model (b) Difference between the full pore-scale solution and the effective solution from the upscaled model Fig. 7 Effective solution from upscaled model and full solution from pore-scale model only two different pore-scale geometries had to be accounted for, but even if more types of different (three-dimensional) geometries are used, the computational gain of rather solving several small, local problems followed by the upscaled problem is significant compared to solving the full pore-scale problem.

Conclusion
In this paper, effective heat transport through thin porous media has been considered. Starting with a pore-scale description, upscaled equations at the Darcy scale have been derived using formal homogenization, while accounting for the porous medium being thin. The Darcy-scale equations rely on the effective permeability as well as the effective heat conductivity of the porous medium which are found through solving local cell problems at the pore scale. These quantities account for the pore geometry as well as the boundary conditions on the top and bottom boundaries of the porous domain.
As a consequence, the effective heat conductivity can be used as an assessment tool for a material's local heat conduction properties without requiring to solve the whole model at the Darcy scale nor at the pore scale. Instead, only small local cell problems need to be solved. Since the cell problems are generally solved on smaller portions of the original domain, they are cheap to solve compared to discretizing and solving the full pore-scale domain. Also, the original pore-scale model is highly coupled, while the cell problems can be solved independently of each other. For homogeneous media, the cell problem solutions and hence effective parameters can be reused, while for heterogeneous media, many cell problems need to be solved. However, solving many small problems is generally cheaper than solving one larger problem from a computational perspective. Furthermore, the cell problems can be straightforwardly solved in parallel. The corresponding upscaled problem is two-dimensional and can be solved on a much coarser grid than the original pore-scale domain. In addition, we have shown that under certain assumptions, the three-dimensional cell problems can be reduced to two-dimensional ones. Hence, the computational costs of solving the cell problems can be further decreased in those cases.
The derived formulation emphasizes the dependence of the effective heat conductivity on the individual heat conductivities of fluid and grain, on the detailed pore geometry as well as on the boundary conditions on the top and bottom boundary. If the grains in a porous medium have a higher heat conductivity than the fluid occupying the void space, we can increase the effective heat conductivity either by decreasing the porosity or, for a constant porosity, by decreasing the minimum distance between the grains, i.e., establishing a better connection of the grains. Furthermore, the application of Dirichlet boundary conditions on the top and bottom boundary has a positive impact on the effective heat conductivities compared to Neumann boundary conditions. This positive impact is in particular strong when the surface area of the more conductive grains toward the top and bottom boundary remains large. However, the effective heat conductivity also depends on the grain shapes and cannot be quantified exclusively by simple parameters such as surface area and porosity.
Hence, to assess the effective heat conductivity of thin porous materials based on their properties at the pore scale, the derived strategy offers equations for determining such effective heat conductivities. These equations are cheap to solve and provide the effective heat conductivity locally, while incorporating the effect from the top and bottom boundary conditions and the pore-scale geometry.