domain: phys spi

text:
M Avila-Lopez

C Robles-Rodriguez

L Tiruta-Barna

A Ahmadi
email: ahmadi@insa-toulouse.fr

Toward thermal autarky for large-scale biogas plants: Dynamic energy modeling for energy efficiency in anaerobic digesters with enhanced multimembrane gasholders

Keywords:

Biogas production via anaerobic digestion could approach thermal autarky (less biogas self-consumption and better energy efficiency) using enhanced solutions for thermal insulation based on new multimembrane gasholder configurations, loss avoidance, and heat recovery strategies. In this study, a predictive and configuration-dependent dynamics energy model is developed to daily and seasonally assess the thermal efficiency and heating requirements of industrial biogas plants involving wet digesters with upper multimembrane gasholders. The model is validated using experimental data from a full-scale biogas plant in operation over a oneyear period.

The energy model involves a set of dynamic energy balances defined for each compartment of the digester and gasholder. All possible heat sources and sinks (ambient air temperature, wind, rain, and solar radiation) and heat exchanges or losses via advection, convection, and conduction, as well as a complete representation of the infrared radiative networks between the surfaces of the digester, are included for each compartment depending on the operating and design parameters. These heat exchanges are subject to fluctuating environmental conditions (e.g., ambient air temperature, wind, rain, and solar radiation).

The results indicate that triple-membrane gasholders with a third insulation membrane made of a suitable material and thickness, together with the involvement of heat recovery from the digestate advective heat, are capable of reducing the overall thermal losses by more than 95 % (e.g., -51 % of the gasholder cover loss and -81 % of the advective digestate loss) and when the waste heat from biogas purification is also valorized, the biogas plant can become thermally self-sufficient.

Introduction

The global need to reduce the use of fossil fuels has triggered the growth of non-fossil bioenergy sources including the production of biogas via anaerobic digestion (AD). Anaerobic digestion has the potential to reduce greenhouse gas emissions by 10 %-13 % and, when used as biomethane, can substitute 26 %-37 % of globally consumed natural gas [START_REF]World Energy Outlook 2018[END_REF]. The benefits of AD include the replacement of conventional chemical fertilizers and the revalorization and management of organic waste [START_REF] Calise | A Review of the State of the Art of Biomethane Production: Recent Advancements and Integration of Renewable Energies[END_REF][START_REF] Holm-Nielsen | The future of anaerobic digestion and biogas utilization[END_REF]. Despite the significant potential of the sector and the existence of more than 100,000 small to large AD units, on average only 2 % of the worldwide potential for AD has been realized [START_REF] Brémond | A vision of European biogas sector development towards 2030: Trends and challenges[END_REF].

Anaerobic digestion is a natural biological process through which different organic materials are degraded into a nutrient-rich digestate and a biogas, which is primarily composed of methane and carbon dioxide. Anaerobic digesters fed on local waste, primarily composed of manure, slurry, waste from food industries, and intermediate crops in co-digestion manure, are categorized as agricultural biogas plants. In France, AD is experiencing significant development, particularly for sewage treatment and for agricultural residues; this is a result of the large potential feedstock deposits in the agricultural sector, representing approximately 70 %-80 % of the potential total feedstock intended for the production of biogas by 2050 [START_REF] Scarlat | A spatial analysis of biogas potential from manure in Europe[END_REF][START_REF] Ademe | La méthanisation en 10 questions[END_REF].

Depending on the content of dry matter (total solids, TS), AD can be performed in wet (<10 %-20 % TS) or dry (>20 %-40 % TS) digesters. Dry AD offers several benefits including low water usage, simple system design, and short residence times. However, a significant number of modern biogas plants are based on wet AD, primarily involving continuously stirred reactors because of their economic performance, linear gas production, and process stability [START_REF] Angelonidi | A comparison of wet and dry anaerobic digestion processes for the treatment of municipal solid waste and food waste[END_REF]. In large-scale wet AD, large stirred digesters operating in the mesophilic regime (~38 • C) are primarily used, with the large digester volumes allowing a maximum expression of the biochemical methane potential of the organic substrate [START_REF] Ademe | Méthanisation. Feuille de route stratégique[END_REF][START_REF] Trame | La méthanisation à la ferme. Guide partique : pour les projets de puissance électrique inférieure à 500 kWe[END_REF].

The produced biogas can be valorized in different ways, generally as electricity and heat production with combined heat and power units, also known as cogeneration, and via biogas purification into biomethane, which is being injected into the gas grid to gradually replace natural gas [START_REF] Holm-Nielsen | The future of anaerobic digestion and biogas utilization[END_REF][START_REF] Scarlat | A spatial analysis of biogas potential from manure in Europe[END_REF][START_REF]Renewables 2021 Global Status Report[END_REF]. Biogas valorization also has the advantage of being energy efficient and environmentally friendly because of its low emission of hazardous pollutants [START_REF] Siciliano | Biogas production from wet olive mill wastes pretreated with hydrogen peroxide in alkaline conditions[END_REF].

Energy efficiency of biogas plants

The overall performance in terms of energy, benefits, and the environmental footprint of an agricultural biogas plant is directly proportional to its energy balance and efficiency, i.e., the amount of energy delivered with respect to the total amount of primary energy produced (the raw biogas from AD), because part of the biogas produced by AD is internally used to heat the digesters [START_REF] Li | Enhancement of methane production in anaerobic digestion process: A review[END_REF][START_REF] Scarlat | A spatial analysis of biogas potential from manure in Europe[END_REF]. This is known as the selfconsumption of the biogas plant, a measure that allows the economic and overall energy performance of a biogas plant to be expressed. Currently, the self-consumption of a biogas plant without thermal insulation or with poor insulation can reach 30 % or even more [START_REF] Bastide | Fiche technique méthanisation[END_REF]; meanwhile, the self-consumption can be from 8 % to 15 % for a digester with insulation, depending on the quality of the digester insulation and the thermal efficiency of the units [START_REF] Ademe | Réaliser une unité de méthanisation à la ferme[END_REF].

The energy efficiency of a biogas plant is presented as the ratio between the final energy produced (electricity and heat from biogas in cogeneration or the production of biomethane) and the primary energy (raw biogas produced in AD). The final energy product of a biogas plant depends on the type of energy targeted, e.g., electricity and heat in cogeneration or biomethane. However, in both cases, the overall energy efficiency of the installation is affected by its thermal efficiency. In cogeneration, the power delivered by the electricity and heat excludes any internal use of biogas in the process. For example, for an AD with a higher thermal requirement, a smaller amount of process heat can be supplied because part of the heat produced is used internally. That is, an increase in the thermal efficiency, and furthermore, the thermal autarky of AD results in an increase in the overall energy efficiency of the biogas plant even in the case of the cogeneration of electricity and heat. Currently, the energy efficiency of AD units with cogeneration is between 50 % and 71 % and the energy efficiency for biomethane production is between 73 % and 85 % [START_REF] Ademe | Suivi technique, économique, et social de 10 installations de méthanisation[END_REF]. The higher energy efficiency values are related to both the electrical efficiency of the cogeneration and the thermal requirements of the AD.

Thermal modeling of large-scale anaerobic digesters and their linked gasholders can play an important role in understanding the conceptual and operational causes of their thermal losses when subjected to fluctuating environmental conditions. Moreover, alternative solutions for thermal insulation, thermal loss avoidance, and heat recovery strategies can be rigorously studied and their potential to approach thermal autarky can be precisely quantified in terms of the self-consumption and energy efficiency of a unit.

Thermal modeling studies of anaerobic digesters vary depending on the depth and specific aims through which the thermal effects and the transfer phenomena are investigated. Approaches range from spatially uniform models [START_REF] Axaopoulos | Simulation and experimental performance of a solar-heated anaerobic digester[END_REF][START_REF] Calise | Modeling of the anaerobic digestion of organic wastes: Integration of heat transfer and biochemical aspects[END_REF][START_REF] Perrigault | Towards thermal design optimization of tubular digesters in cold climates: A heat transfer model[END_REF] to the use of computational fluid dynamics (CFD) to consider spatial heterogeneities [START_REF] Liu | Investigation on the heat loss characteristic of underground household biogas digester using dynamic simulations and experiments[END_REF][START_REF] Rynkowski | Heat loss analysis in the semi-buried anaerobic digester in Northeast Poland[END_REF][START_REF] Wu | Development of 3-d anaerobic digester heat transfer model for cold weather applications[END_REF]. Spatially uniform models are better adapted to analysis and optimization studies over long dynamic periods, especially in large-scale industrial digesters, because the CFD approach is computationally expensive and very complex when applied to a non-Newtonian digestate in the presence of bubbles and suspended solids [START_REF] Hreiz | Modeling and simulation of heat transfer phenomena in a semi-buried anaerobic digester[END_REF].

There are several compartmental approaches in which spatially

Nomenclature

A Surface area (m 2 ) A 0 Amplitude of the surface-temperature fluctuation (K) c p Heat capacity (J kg

(J kg -1 ) h Convective heat transfer coefficient (W m -2 K -1 ) h ′ Specific enthalpy (J kg -1 ) h Height of an equivalent cuboid (m) J Radiosity (W m -2 ) K T Insolation clearness index (-) k Thermal conductivity (W m -1 K -1 ) L c
Characteristic length (m) ṁ

Mass flowrate (kg s -1 ) n Day of the year (-)

P Perimeter (m) R Thermal resistance (K W -1 )

S eq

Base surface of an equivalent cuboid (m ρ

U Thermal transmittance (W m -2 K -1 ) u Flow speed (m s -1 ) V Volume (m -3 ) ẆM Dissipated heat
Density (kg m -3 ) σ Stefan-Boltzmann constant (W m -2 K -4 ) ψ Ground albedo (-)
uniform heat models have been used to quantify the thermal losses of AD given different digester configurations, such as rural tubular digesters [START_REF] Perrigault | Towards thermal design optimization of tubular digesters in cold climates: A heat transfer model[END_REF] and industrial cylindrical digesters [START_REF] Hreiz | Modeling and simulation of heat transfer phenomena in a semi-buried anaerobic digester[END_REF][START_REF] Calise | Modeling of the anaerobic digestion of organic wastes: Integration of heat transfer and biochemical aspects[END_REF][START_REF] Vilms Pedersen | Management and design of biogas digesters: A non-calibrated heat transfer model[END_REF]. Hreiz et al. [START_REF] Hreiz | Modeling and simulation of heat transfer phenomena in a semi-buried anaerobic digester[END_REF] provided a detailed assessment of the heat losses of a medium-scale semi-buried agricultural digester with a single-membrane gasholder using zero-dimensional thermal modeling as a function of the operating and weather conditions; in this study, attempts were made to better represent the solar heating and to predict the biological heat of digestion according to the composition of the input organic matter. Calise et al. [START_REF] Calise | Modeling of the anaerobic digestion of organic wastes: Integration of heat transfer and biochemical aspects[END_REF] coupled thermal and biochemical models to examine a doublemembrane AD plant via a simplified ADM1 model for the biological processes and the digester heat balance for the heat modeling; in this case, constant convective heat transfer coefficients were used for the losses from the biogas and the air gap while the losses from the cover and the inner membrane were neglected. Vilms Pedersen et al. [START_REF] Vilms Pedersen | Management and design of biogas digesters: A non-calibrated heat transfer model[END_REF] developed a thermal model for a fixed-dome digester and an industrial anaerobic digester using a thermal resistance network approach that calculates the digestate temperature according to only the ambient and operating conditions. Available heat models usually involve good representations of the heat transfer via advection, conduction, and convection for singlemembrane gasholders or gasholders modeled as a single-layer heat barrier but lack (i) a detailed geometric representation of the gasholder layers and their complete energy balance and associated heat transfer, (ii) a complete internal longwave infrared radiative heat exchange network between the gasholder membranes, and (iii) a dynamic prediction of the convective and radiative losses in the gasholder according to the geometric considerations and the varying operating conditions.

Goals and novelties

The first goal of this study is to develop and experimentally validate a predictive tool to dynamically evaluate the thermal efficiency and heating requirements of large-scale anaerobic digesters with multimembrane gasholders to provide a reliable thermal assessment of a digester over a long operational period. The model will serve as a predictive tool for sensitivity analyses to optimize the design and operating conditions of anaerobic digesters. The second goal is to evaluate the potential of new large-scale anaerobic digesters to approach thermal autarky when multimembrane gasholders and various heat recovery strategies are applied. The main goal in this respect is to propose and justify viable solutions via energy modeling to enable modern largescale digesters to reduce their thermal losses.

Accordingly, the primary and distinctive characteristics of this study and the developed dynamic energy model are as follows. (1) Experimental data are collected from a large-scale biogas plant over an entire year to validate the energy model and to identify feasible solutions allowing thermal autarky to be approached by a real-world industrial biogas plant. (2) The digester and multimembrane gasholder are divided into multiple internal compartments including the digestate, walls, biogas, membranes, gasholder cover, and air ventilation zone, depending on the number of gasholder membranes. (3) The modeling is flexible and configuration dependent, capable of considering gasholders with single-to triple-membrane configurations, and varying membranes and geometries can be applied with adjustable thermo-optimal properties. (4) The advective, conductive, convective, and radiative heat exchanges for each compartment are fully applied together with all dynamic environmental conditions, such as the air and soil temperatures, solar radiation, wind, and rain. [START_REF] Scarlat | A spatial analysis of biogas potential from manure in Europe[END_REF] The longwave infrared radiative heat losses are fully included via a radiative network approach defining the surface radiosities of each modeling compartment. (6) Dynamic soil temperature modeling is used according to the soil depth and the division of the soil heat losses into the aboveground and underground wall losses to model the effect of buried and semi-buried digesters on the thermal efficiency. [START_REF] Angelonidi | A comparison of wet and dry anaerobic digestion processes for the treatment of municipal solid waste and food waste[END_REF] The global solar irradiation with clouds is dynamically estimated using the geographically correlated sky clearness index to determine the net reliable shortwave solar radiation reaching the digester. ( 8) All design and operating parameters for the energy model are exhaustively incorporated to form a powerful predictive tool for the thermal design and optimization of AD.

Method -dynamic energy modeling of large-scale AD units

The general design of modern semi-buried wet anaerobic digesters at farms equipped with multimembrane gasholders is illustrated in Fig. 1. The top membrane acts as a protective cover, and below this cover, supplementary membranes may be included, either for storage of the produced biogas or for thermal insulation. As shown in Fig. 1, in modern biogas plants, the multimembrane gasholder is placed above the digester. The biogas produced in the digester is stored in the overhead gas above the digestate by means of a flexible storage membrane (inner membrane), which is then separated from the surrounding air by an air pocket and a protective membrane covering the roof (Fig. 1). To adjust the inner membrane according to the need to store or destock the biogas, the contained air is increased via an air blower or is expelled via a damper valve. Insulating rubber with a mechanical consistency such as polyvinyl chloride or ethylene propylene diene monomer rubber is generally used with a thickness greater than 1.5 mm. The role of the inner flexible membrane associated with the air blower is to maintain a relative biogas pressure between 15 mbar and 30 mbar.

In the case of industrial wet digestion under mesophilic conditions, digesters are usually well-mixed continuous stirred tank reactors. To consider heat exchanges with the environment, the weather conditions need to be adapted to geographically suit the actual location of the digester. The main material inputs to the system are influents composed of agricultural wastes, slurry, and liquid digestate, while the main material outputs are the digestate and the produced brute biogas. The main energy inputs are the heat used to maintain the mesophilic temperature inside the digester and the electricity used for the stirrer.

To consider the heat exchanges between the different components of digesters with multimembrane gasholders, the unit is divided into several compartments, the temperatures of which are considered as being uniform throughout the compartment but variable over time. For example, in the case of a double-membrane gasholder, the compartments are the (i) digestate (d), (ii) biogas (b), (iii) air ventilation zone (air), (iv) inner membrane (im), (v) cover (c), and (vi) walls (w). A double-membrane gasholder (with an inner membrane and a cover) is adjusted using air pressure via a blower and a pressure regulator to maintain the shape and structure of the cover over the gasholder.

If the unit is equipped with two insulation membranes above the gasholder, a total of six energy balances per digester can be written (Eqs.

(1)-( 6)), one energy balance per compartment. The energy balance equations are derived from the advective heat transfers and the appropriate set of conductive, convective, and radiative heat transfers. The radiative transfers include heat transfer from the digestate to the cover ( Qd,net(rad) ), radiation from the cover to the atmosphere ( Qc,at(rad) ), and solar irradiance absorbed by the cover ( QSunc(rad) ), as well as all infrared radiative exchanges between surfaces, where the geometries and dimensions are considered via the relative view factors. In Section 2.2, all other terms related to the reciprocal heat transfer phenomena between the compartments are completely defined using the relevant theoretical and empirical correlations.

According to the heat balance equations, the main dynamic variables of the system per digester are the temperatures (T d : digestate, T b : biogas, T w : wall, T im : inner membrane, T air : air ventilation zone, and T c : cover) when the energy supplied by the heat exchanger to heat the digester is known ( QHE ). Note that, when dynamic experimental data concerning the temperature in the digestate (T d ) are available via measurements, Eq. ( 1) can easily be modified such that the heating energy supplied by the heat exchanger ( QHE ) becomes a new variable that can be predicted over time.

Digestate energy balance (d)

ρ d V d c p,d dT d dt = Δ Ḣd -Δ Ḣd,b -Qd,b -Qd,w -Qd,net(rad) + ẆM -Qd,gr + QHE -Δ Ḣr (1)
Biogas energy balance (b)

ρ b V b c p,b dT b dt = Δ Ḣd,b -Δ Ḣb -Qb,im + Qd,b -Qb,w (2)
Air ventilation zone energy balance (air)

ρ air V air c p,air dT air dt = Δ Ḣair + Qim,air -Qair,c (3)
Inner membrane energy balance (im)

ρ im V im c p,im dT im dt = Qb,im -Qim,air + Qim,inf,net(rad) -Qim,sup,net(rad) (4)
Cover energy balance (c)

ρ c V c C c
dT c dt = Qair,c -Qc,at -Qc,rain + Qc,inf,net(rad) -Qc,at(rad) + QSunc(rad) [START_REF] Scarlat | A spatial analysis of biogas potential from manure in Europe[END_REF] Wall energy balance (w) ρ w V w c p,w dT w dt = Qd,w + Qb,w -Qw,at + Qw,net(rad) -Qw,at(rad) + QSunw(rad) [START_REF] Ademe | La méthanisation en 10 questions[END_REF] The same modeling strategy can be applied to units with digesters equipped with a single-or triple-membrane gasholder. In the former case, Eqs. ( 3) and (4) can be removed to allow the biogas compartment to be directly in contact with the cover of the gasholder, resulting in a set of four energy balance equations (i.e., for the cover, biogas, digestate, and walls). In the latter case, two new energy balances need to be added, one representing the new insulation membrane (similar to Eq. ( 4)) and the other to balance the energy in the new gas compartment appearing between the inner membranes (similar to Eq. ( 2)), resulting in a set of eight energy balance equations (i.e., for the cover, membrane 1, membrane 2, biogas 1, biogas 2, ventilation air, digestate, and walls).

Heat transported through advective flows

The advective heat flows, designated by Δ Ḣ in the energy balance equations (Eqs. (1)-( 6)), consist of the amount of energy transported as inputs or outputs by the mass flow in the digester and only depend on the enthalpy and the flowrate [START_REF] Smith | Introduction to chemical engineering thermodynamics[END_REF]. According to the energy balance equations, these advective flows are the overall enthalpy flow between the digester feed and the outgoing digestate (Δ Ḣd ), the enthalpy flow of biogas leaving the digestate (Δ Ḣd,b ), the enthalpy flow of biogas leaving the gasholder (Δ Ḣb ), and the overall enthalpy flow of air ventilated below the gasholder cover (Δ Ḣair ).

Δ Ḣd = Δ Ḣf -Δ Ḣd,out = ṁf h ′ f,T f -ṁd h ′ d,T d (7) Δ Ḣd,b = ṁb h ′ b,Td (8)
Δ Ḣb = ṁb h ′ b,T b with h ′ b = h ′ biogas + x H2O Δh vap,H2O (9)
Δ Ḣair = Δ Ḣair,in -Δ Ḣair,out = ṁair ( h ′ air,in -h ′ air,out ) (10)
Here, Δ Ḣf and Δ Ḣd,out are the enthalpy flows of the incoming influent and outgoing liquid digestate, respectively, ṁ indicates the mass flowrate, and h ′ indicates the specific enthalpy determined according to the temperatures and compositions of the components and mixtures.

The thermodynamic properties of the biogas are evaluated dynamically given the biogas temperature and composition, where the biogas is assumed to be an ideal incompressible gas. Knowing the temperature, the water content in the biogas is determined by calculating the saturated vapor pressure of water. The thermodynamic properties of the digester influents, as well as the liquid digestate, are primarily composed of water and are therefore treated the same as liquid water (because of their low dry matter composition) and vary according to the temperature. The air is also assumed to be an ideal gas with temperaturedependent thermodynamic properties, entering the ventilation zone at ambient temperature (T amb ) and exiting this zone at a higher temperature.

The thermodynamic properties (i.e., the viscosity, density, specific heat, and thermal conductivity) of the digestate, biogas (CH 4 , CO 2 , and H 2 O) and air, as well as the saturated vapor and liquid enthalpy of water in the biogas, are temperature dependent and calculated at an operating pressure using the CoolProp library [START_REF] Bell | [END_REF]. The thermal expansion coefficient of the digestate is calculated using the freshwater library [START_REF] Robb | PyPI[END_REF]. Note that, throughout the dynamic energy modeling, all temperaturedependent properties, such as the thermal conductivities and saturation properties, as well as the convective heat transfer coefficients are dynamically re-evaluated based on the updated values of the state variables.

Convective and conductive heat transfers 2.2.1. Convective heat transfer from digestate to biogas

The heat flow via convection from the digestate to the biogas is calculated as follows:

Qd,b = h d,b A d (T d -T b ) (11)
where h d,b is the convective heat transfer coefficient between the digestate and the biogas and A d is the digestate free surface area. The biogas flow inside the digester and between the digestate-biogas interface is assumed to be very low, such that the convection above the digestate is assumed to be natural. Therefore, h d,b is determined following the principle of convection on the upper side of a horizontal hot plate using Eqs. (A.6) and (A.

Convective and conductive heat transfer from digestate to aboveground walls

The heat flow from the digestate through the aboveground walls is determined by

Qd,w = A w,overground (T d -T w ) 1 hd,w + ∑ n i=1 δi ki ( 12
)
where h d,w is the convective heat transfer coefficient between the digestate and the walls, A w,overground is the aboveground wall surface area in contact with the digestate, δ i is the thickness of the walls and insu-lation, and k i is the thermal conductivity of the walls and wall insulation.

The flow movement of the digestate is assumed to be slow; therefore, h d,w is treated as natural convection and is determined using the Churchill and Chu [START_REF] Churchill | Correlating equations for laminar and turbulent free convection from a vertical plate[END_REF] correlation for natural convection in a vertical plate (Eqs. (A.9) and (A.10)).

Conductive heat transfer from digestate to underground walls and ground

The soil temperature is influenced by many different factors and changes in response to the dynamic meteorological conditions, soil properties, geographic location, and human management. All of these factors make it difficult to define the thermal regime of the soil properties.

Different correlations and methods to estimate the heat flux from the digestate through the underground walls and the bottom of the digester to the ground can be found in the literature; these include thermal resistance approaches [START_REF] Perrigault | Towards thermal design optimization of tubular digesters in cold climates: A heat transfer model[END_REF], as well as transient semi-infinite solid [START_REF] Gebremedhin | Heat transfer model for plugflow anaerobic digesters[END_REF] and two-dimensional transient CFD [START_REF] Hreiz | Modeling and simulation of heat transfer phenomena in a semi-buried anaerobic digester[END_REF] simulations. This study uses the EN ISO 13,370 method for heat transfer for building elements in contact with the ground [START_REF]Thermal performance of buildings -Heat transfer via the ground -Calculation methods[END_REF]. This method has been standardized using reliable real-world representations of thermal resistances and transfers to walls and the ground and allows the temperature gradient of the soil with depth to be considered, as well as the dynamic influence of the ambient temperature and the surface temperature of the soil in the upper soil layer.

The heat flux from the digestate to the ground is calculated as follows:

Qd,gr = ( A floor U floor (T d -T s (z, t) ) ) + (A wu U wu (T d -T s (z, t) ) ) (13)
where A floor is the bottom surface area of the digester, A wu is the wall surface area in contact with the ground, U floor is the thermal transmittance of the bottom floor, U wu is the thermal transmittance of the underground walls, T s (z, t) is the soil dynamic temperature at the depth of the digester, and T s,z is the mean soil temperature between the soil surface and the depth of the digester. U floor and U wu are calculated using Eqs. (A. [START_REF] Axaopoulos | Simulation and experimental performance of a solar-heated anaerobic digester[END_REF]) and (A.17 Because the soil temperature (T s ) is a function of time (t) and depth (z) in the soil, it oscillates with a sinusoidal shape over time, with the amplitude of its variation being determined according to the minimum and maximum ambient air temperatures over the chosen time period [START_REF] Hillel | Introduction to Environmental Soil Physics[END_REF]:

T s (z, t) = T s,z0 + T amb,max -T amb,min 2 • exp ( -z ̅̅̅̅̅ ̅ π dα √ ) • cos ( 2π d [t -t 0 ] -z ̅̅̅̅̅ ̅ π dα √ ) (14)
Here, T s,z0 is the mean soil temperature at the surface; T amb,max and T amb,min are the maximum and minimum ambient air temperatures, respectively, during the chosen time period, t; z is the depth; d is the duration of the chosen period (e.g., in days or hours); α is the thermal diffusivity of the soil; and t 0 is the time delay from an arbitrary starting time of the period (generally the first day, hour, or second of the period) to the appearance of the minimum temperature in the chosen period. In Eq. ( 13), T s (z, t) is set to the temperature at the bottom depth and T s (z, t) is set to the average soil temperature between the surface T s (0, t) and the digester base T s (z, t).

Convective and conductive heat transfer from biogas to walls

Equation ( 15) gives an expression of the global heat flow from the biogas to the aboveground walls by merging the conduction through the walls and the convection in the biogas according to their respective thermal resistances.

Qb,w =

A w,over-bg (T b -T w )

1 hb,w + ∑ n i=1 δi ki (15)
Here, h b,w is the convective heat transfer coefficient between the biogas and the walls and A w,over-bg is the aboveground wall surface area in contact with the biogas. Similar to the heat transfer between the digestate and the walls, natural convection is assumed and h b,w is obtained using Eqs. (A.9) and (A.10) for natural convection in a vertical plate [START_REF] Churchill | Correlating equations for laminar and turbulent free convection from a vertical plate[END_REF].

Convective heat transfer from walls to ambient air

The heat flow from the aboveground walls to the ambient air is formulated as

Qw,at = h w,at A wo (T w -T amb ) (16)
where h w,at is the convective heat transfer coefficient between the walls and the ambient air and A wo is the total aboveground wall surface area in contact with the ambient air. Depending on the wind speed and its flow regime, the heat loss to the air is either modeled as natural convection or as forced convection. In this study, two relevant situations are distinguished: (1) no wind, where the losses are dominated by natural convection and Eqs. (A.9) and (A.10) are used, and (2) the presence of wind, where the Churchill and Bernstein [START_REF] Churchill | A correlation equation for forced convection from gases and liquids to a circular cylinder in crossflow[END_REF] correlation (Eq. (A.13)) for the convection of air flow normal to a circular cylinder in a cross flow is used such that the entire Reynolds and Prandtl number ranges are covered.

Convective and conductive heat transfer from biogas to inner membrane

The heat exchange between the biogas and the inner membrane is formulated in Eq. ( 17) as the sum of the convective and conductive heat flows.

Qb,im =

A im (T b -T im )

1 hb,im + ∑ n i=1 δi ki (17)
Here, h b,im is the convective heat transfer coefficient between the biogas and the inner membrane and A im is the inner membrane surface area. The convection is assumed to be natural, and the Nusselt number is determined using the correlation provided in Appendix A.2.1 for horizontal plates [START_REF] Mcadams | Heat Transmission[END_REF]: Eqs. (A.6) and (A.7) for T b > T im represent convection on the lower side of a horizontal cold plate and Eq. (A.8) for T b ≤ T im describes convection on the lower side of a horizontal hot plate.

Convective heat transfer from inner membrane to air ventilation zone

The heat flow from the inner membrane to the air ventilation zone is represented by

Qim,air = h im,air A im (T im -T a ) ( 18
)
where h im,air is the convective heat transfer coefficient between the inner membrane and the airflow. Depending on the air speed and the air flow, the convection is assumed to be forced, natural, or mixed. The mean air speed (u air ) is estimated from the volumetric flowrate of the air ( V) and the mean cross-sectional area of the air ventilation zone. As illustrated in Fig. 2, assuming a flow inside a cuboid with a square base with a surface area of S eq (equal to the real surface of the membrane), the mean rectangular cross-sectional area between the inner membrane and the cover Fig. 2. Theoretical configuration for the determination of the air speed.

can be obtained as the ratio between the actual volume of the air ventilation zone and the actual surface area of the membrane (h = V airgap /S eq ). Therefore, the mean cross section is V airgap / ̅̅̅̅̅̅ ̅ S eq √ ; this parameter is required for the approximation of the mean air speed across the ventilation zone.

If Gr/Re 2 ≥ 10, where Gr is the Grashof number and Re is the Reynolds number, natural convection takes place and the Nusselt number is calculated using Eqs. (A.6) and (A.7) for T im > T air and Eq. (A.8) for T im < T a . If Gr/Re 2 ≤ 0.1, the heat transfer is described by forced convection and the forced convection correlation in a parallel flow (Eqs. (A.11) and (A.12)) is used. If 0.1 < Gr/Re 2 < 10, mixed convection occurs and the Nusselt number is obtained as follows [START_REF]Heat Exchanger Desing Handbook[END_REF]:

Nu n Mix = Nu n N + Nu n F (19)
Here, Nu N and Nu F refer to the Nusselt numbers for natural and forced convection, respectively. The value of the coefficient n depends on the geometry and the flow direction; for transverse flows involving horizontal plates, a value of n = 3.5 is recommended [START_REF] Incropera | Fundamentals of heat and mass transfer[END_REF].

Convective and conductive heat transfer from air ventilation zone to cover

The global heat flow from the ventilation air compartment to the cover is calculated using the following equation:

Qair,c = A c (T air -T c ) 1 hair,c + ∑ n i=1 δi ki ( 20
)
where h air,c is the convective heat transfer coefficient between the air gap and the cover and A c is the cover area. Similar to the heat transfer between the inner membrane and the pressurized air, the convection can be natural, forced, or mixed depending on the value of Gr/Re 2 . The same criteria as given in Section 2.2.7 are used to determine if natural, forced, or mixed convection is applied, using Eqs. (A.6) and (A.7) for T a > T c and Eq. (A.8) for T a < T c under natural convection conditions. In the event of forced or mixed convection, Eqs. (A.11) and (A.12) and Eq. ( 19) are used, respectively.

Convective heat from cover to ambient air

Equation ( 21) describes the heat transfer from the cover to the ambient air via convection.

Qc,at = h c,at A c (T c -T amb ) (21)
The convective heat transfer coefficient in this equation (h c,at ) depends on the wind speed (u wind ). If there is no wind (u wind = 0), natural convection is assumed where the corresponding Nusselt numbers are determined using Eqs. (A.6) and (A.7) when T c > T amb and Eq. (A.8) when T c < T amb . Conversely, in the presence of wind (u wind ∕ = 0), the forced convection correlations are applied as given in Eqs. (A.11) and (A.12).

Convective heat loss via rain

The additional heat loss resulting from the flow of rain over the cover is determined by introducing a new convective heat loss coefficient (h c,rain ) reflecting the appearance of a thin water layer over the external cover area as given in Eq. [START_REF] Hreiz | Modeling and simulation of heat transfer phenomena in a semi-buried anaerobic digester[END_REF].

Qc,rain = h c,rain A c (T c -T amb ) (22)
In the presence of rain, forced convection is considered using Eqs. (A.11) and (A.12). The volumetric rain flowrate over the cover is estimated from the total precipitation data and the horizontal projected surface of the cover. A liquid film of water with a thickness of 1 mm on the cover is assumed to represent the presence of a thin water layer on the cover where the dripping speed of the rain and the total amount of rainfall are averaged on a monthly basis.

Radiative heat transfer 2.3.1. Radiation exchange between internal surfaces

The different longwave infrared heat fluxes inside the confined volumes of the digester and gasholder, i.e., the radiative heat exchanges between the digestate, walls, and membranes from one side and those between the insulation and protective membranes are precisely addressed in the radiation network approach via the surface radiosity variables (J i ) [START_REF] Oppenheim | Radiation Analysis by the Network Method[END_REF][START_REF] Rohsenow | Handbook of Heat Transfer[END_REF]. The net radiative flow leaving any surface in a confined space can be written as follows:

q i = ε i σT 4 i -α i J i (1 -α i )/A i = ∑ n j=1 ( A i F ij )( J i -J j ) ( 23
)
where σ is the Stefan-Boltzmann constant, ε i is the emissivity, α i is the absorptivity, T i is the temperature, and A i is the area of the surface i. J i and J j are the radiosities from the surfaces i and j, respectively, and F ij is the view factor between the surfaces i and j. Equation ( 23) represents the net radiative flow leaving the surface/ node i considering other surfaces/nodes (j) in a given confined space and receiving their respective radiosities (J j ). In this equation, the net radiative flow leaving any surface (i) is the difference between the total radiative emissions according to the surface temperature (T i ), emissivity (ε i ), and radiation absorbed by the surface based on its absorptivity (α i ); the absorbed radiation is the total incoming radiation from the other surfaces minus the reflected radiation from the surface (i). A demonstration of how Eq. ( 23) can be obtained as a result of the balance between the absorptive, emissive, and reflective radiation over any surface in a confined volume is given in Supplementary material (Section 1). The calculation of the view factors between the surfaces, F ij , is presented in Supplementary material (Section 2).

According to the radiation network approach as formulated in Eq. ( 23), the net radiative heat flow leaving any surface inside the digester and gasholder can be formulated as shown in Eqs. ( 24)- [START_REF] Churchill | Correlating equations for laminar and turbulent free convection from a vertical plate[END_REF].

Net radiative flow from the digestate surface

Qd,net(rad) = ε d σT 4 d -α d J d (1 -α d )/A d = ( A d F d,im,inf )( J d -J im,inf ) + ( A d F d,w ) (J d -J w ) (24)
Here, Qd,net(rad) is the net radiation leaving the digestate surface when viewing the inferior side of the inner membrane (inf) and the walls given their radiosities (J im,inf and J w , respectively) and ε d and α d are the emissivity and absorptivity of the digestate, respectively, which are equal according to Kirchoff's law and are set to 0.9 in agreement with the literature [START_REF] Kishore | A heat-transfer analysis of fixed-dome biogas plants[END_REF][START_REF] Rennie | Liquid manure storage temperature is affected by storage design and management practices -A modelling assessment[END_REF].

Net radiative flow from the inferior side of the inner membrane

Qim,inf,net(rad) = ε im,inf σT 4 im -α im,inf J im,inf ( 1 -α im,inf )/ A im = ( A im F im,inf ,d )( J im,inf -J d ) + ( A im F im,inf ,w )( J im,inf -J w ) (25)
Here, Qim,inf,net(rad) is the net radiation leaving the inferior surface of the inner membrane when viewing the digestate and the walls given their radiosities (J d and J w , respectively) and ε im,inf and α im,inf are the emissivity and absorptivity, respectively, of the inferior surface of the inner membrane corresponding to the membrane material properties.

Net radiative flow from the interior side of the walls

Qw,net(rad) = ε w,int σT 4 w,int -α w,int J w ( 1 -α w,int )/ A w = ( A w F w,d ) (J w -J d ) + ( A w F w,im,inf )( J w -J im,inf ) (26)
Here, Qw,net(rad) is the net radiation leaving the wall surfaces when viewing the inferior side of the inner membrane and the digestate given their radiosities (J im,inf and J d , respectively) and ε w,int and α w,int are the emissivity and absorptivity, respectively, of the interior side of the walls corresponding to the optical properties of the concrete. T w,int represents the wall temperature inside the digester, which is in contact with the biogas.

Net radiative flow from the superior side of the inner membrane

Qim,sup,net(rad) = ε im,sup σT 4 im -α im,sup J im,sup ( 1 -α im,sup )/ A im = ( A im F im,sup,c )( J im,sup -J c ) (27)
Here, Qim,sup,net(rad) is the net radiation leaving the superior side of the inner membrane (sup) when viewing the inferior side of the cover given its radiosity (J c ) and ε im,sup and α im,sup are the emissivity and absorptivity, respectively, of the superior surface of the inner membrane.

Net radiative flow from the inferior side of the cover

Qc,net(rad) = ε c,inf σT 4 c -α c,inf J c ( 1 -α c,inf )/ A c = ( A c F c,im,sup )( J c -J im,sup ) (28)
Here, Qc,net(rad) is the net radiation leaving the inferior part of the cover when viewing the superior side of the inner membrane given its radiosity (J im,sup ) and ε c,inf and α c,inf are the emissivity and absorptivity, respectively, on the inferior side of the cover.

Net longwave infrared radiation between walls and surroundings

The longwave infrared radiation between the walls and the surroundings is calculated as the net radiation exchanges from the lateral aboveground walls on one side and the ground, sky, and surroundings considered as radiative bodies at ambient temperature on the other side [START_REF] Mcadams | Heat Transmission[END_REF][START_REF] Walton | Thermal Analysis Research Program Reference Manual[END_REF]:

Qw,at(rad) = F w-sky ε w,ext σA w ( T w 4 -T Sky 4 ) + F w-gr ε w,ext σA w ( T w 4 -T amb 4 ) + F w-sur ε w,ext σA w ( T w 4 -T amb 4 ) (29)
where ε w,ext is the emissivity related to the insulation material used on the exterior side of the walls and T Sky is the effective sky temperature, which is lower than the ambient air temperature determined via Eq. [START_REF] Hillel | Introduction to Environmental Soil Physics[END_REF]. The view factors from the walls to the sky ( F w-sky ) , from the walls to the surroundings at ambient air temperature (F w-sur ), and from the walls to the ground ( F w-gr ) are illustrated in Fig.

3 [START_REF] Walton | Thermal Analysis Research Program Reference Manual[END_REF]. Different correlations can be found in the literature to determine the sky temperature (T Sky ) for models based on clear to cloudy sky condi-tions. Regarding the energy model presented in this study, it is essential to choose a sky temperature model involving the cloudiness effect to determine the sky temperature because such models are more relevant to European sky conditions. Clouds usually increase the effective sky temperature, causing it to approach the ambient air temperature [START_REF] Aigarni | Survey of sky effective temperature models applicable to building envelope radiant heat transfer[END_REF]; this highlights the significant effect of cloudy weather on T Sky . Here, the empirical correlation proposed by Fuentes [START_REF] Fuentes | A simplified thermal model for flat plate photovoltaic arrays[END_REF] for cloudiness as a function of the ambient temperature was chosen:

T Sky = 0.037536T 1.5 amb + 0.32T amb (30)

Net longwave infrared radiation between cover and sky

The net radiative heat flow between the cover and the sky is formulated in Eq. [START_REF] Hillel | Introduction to Environmental Soil Physics[END_REF] assuming that the sky is a gray body emitting longwave infrared radiation according to its temperature T Sky .

Qc,at(rad) = ε co,sup σA c ( T c 4 -T Sky 4 ) (31)
Here, ε co,sup is the emissivity of the superior side of the cover; this value depends on the optical properties of the materials used for the protective cover.

Shortwave solar irradiance absorbed by the cover

The shortwave solar irradiance absorbed by the cover, as formulated in Eq. ( 32), depends on the instantaneous global solar irradiation on a surface at the Earth's surface (G i H ), the solar absorptivity of the surface materials (α s,co,sup ), and the cover area (A c ). The component that is received from the sky is proportional to the view factor between the cover and the sky (F c-sky ), and the component resulting from the reflection of the ground is determined by the corresponding view factor between the cover and the ground (F c-gr ), as well as the soil albedo (ψ).

QSun-c(rad) = G i H α s,co,sup A c ( F c-sky + F c-gr ψ ) (32)
Here, the view factors between the cover and the sky (F c-sky ) and between the cover and the ground (F c-gr ) are determined using the method described in Section 2.3.2 but with θ set to the tilt angle of the cover.

As described in Supplementary material (Section 3), the global solar irradiation on the ground (G H ) is the sum of two shortwave components: the beam irradiation (B H ) and the diffuse radiation (D H ). In Supplementary material (Section 3), G H is first corrected to include the effect of cloudiness (using the daily clearness index regression K T , which is specific to the location of the plant) and is then converted to the instantaneous G i H using the solar hour angle (ω) for the location of the plant and its respective geographic coordinates.

Shortwave solar irradiance absorbed by the walls

The shortwave solar irradiance absorbed by the walls is also determined based on the global solar irradiation on a surface at the Earth's surface (G i H ) but with view factors, optical properties, and surfaces relative to the lateral walls:

QSun-w(rad) = G i H α s,w,ext A w,over ( F w-sky + F w-gr ψ ) ( 33
)
where α s,w,ext is the solar absorptivity of the exterior side of the walls, which is set to 0.74 corresponding to dark green color properties. ψ corresponds to the ground albedo and has values varying from 0.13 to 0.16 depending on the month. The term G i H times F w-sky represents the contribution of the direct solar irradiance and the diffuse irradiance, and the term G i H times F w-gr ψ represents the contribution of the reflected ground irradiance.

Other heat terms 2.4.1. Heat of biodegradation

The heat of biodegradation (Δ Ḣr ) is estimated using the approach of Lindorfer et al. [START_REF] Lindorfer | Self-heating of anaerobic digesters using energy crops[END_REF], assuming that the organic substrates can be represented as a mixture of three main component categories: (i) carbohydrates represented by glucose (C 6 H 12 O 6 ), (ii) fats represented by palmitic acid (C 16 H 32 O 2 ), and (iii) proteins represented by alanine (C 3 H 7 NO 2 ). Using this representation, average digester feed (solid waste plus cattle slurry feed) is composed of 2.26 wt% lipids, 14.91 wt% carbohydrates, and 4.81 wt% proteins. The heat of biodegradation is determined by incorporating the corresponding enthalpies of biodegradation for carbohydrates (glucose), fats (palmitic acid), and proteins (alanine), these values being -768.7 kJ kg -1 , 2123.4 kJ kg -1 , and 2228.1 kJ kg -1 , respectively. The global heat of biodegradation can then be determined accordingly to the mass fraction of each subcategory of organic matter (w i ), the respective specific biodegradation enthalpy (Δ Ḣi ), and the total mass flowrate of the solid waste plus the slurry ( ṁsw ).

Δ Ḣr = ṁsw ∑ w i Δ Ḣi (34)

Integration and solution procedure

Despite the large number of properties and parameters and their respective equations, the main system of equations for a configuration composed of a digester and a double-membrane gasholder can be summarized into six differential balance equations per digester (Eqs. ( 1)-( 6)) and five algebraic equations describing the infrared radiative network between the internal surfaces in each digester (Eqs. ( 24)-( 28)). The global mathematical problem is therefore a system of differential algebraic equations involving a total of 11 equations and 11 independent variables per digester, i.e., the 6 temperatures given in Section 2 and the 5 internal surface radiosities (J d , J c , J w , J im,sup , and J im,inf ).

The differential equations were integrated with a variable step size in Python 3.7 via the built-in ordinary differential equation solver function scipy.integrate.ode using the lsoda integrator solver with a method based on backward differentiation formulas [START_REF]The SciPy community[END_REF]. The scipy.optimize. fsolve library was used in parallel to resolve the algebraic equations [START_REF]The SciPy community[END_REF].

Because the dynamic variation of the digestate temperature is a commonly monitored parameter that is recorded in all industrial biogas plants to regulate the AD temperature, the digestate temperature is considered to be a known dynamic parameter in this study. Therefore, the heat supply required from the heat exchanger ( QHE ) is considered to be a system variable. The QHE value calculated from the model represents the predicted heat demand of a given digester; this enables the further analysis of the energy performance of industrial biogas plants.

Results and discussions

Case study: Application to large-scale on-farm anaerobic digesters with multimembrane gasholders

A semi-buried wet anaerobic digester and post-digester in series with double-membrane gasholders, referred to here as Biometharn, located in France and intended for biogas injection, was used to assess the thermal model and its experimental validation. The digester and post-digester are presented graphically in Fig. 4. The digesters are industrial-scale well-mixed continuous stirred tank reactors. The digester and postdigester have walls and bottoms with diameters of 25 m and 18 m, respectively. The double-membrane gasholder has an inner membrane made of an ethylene propylene diene monomer and an outer membrane or cover made of polyvinyl chloride. The digester and post-digester are buried 2.5 m belowground. The cover protects the inner membrane and maintains the structure, and the inner membrane stores the biogas. The digestate operating temperature is slightly variable but is globally stable close to 40 • C. The digester is equipped with two three-blade propeller mixers and one four-blade paddle agitator. The post-digester is equipped with a single three-blade propeller mixer. The characteristics of both the digester and the post-digester (e.g., dimensions, shapes, and materials) are listed in Table 1. The mechanical mixing system for the digester is powered with a 15-kW motor, and that for the post-digester is powered with a 4-kW motor. Considering an efficiency of 90 %, the dissipated heat supplied by the work of the mixing systems ( ẆM ) was assumed to be 13.5 kW and 3.6 kW for the digester and post-digester, respectively (Table 1).

The walls are made of 20-cm-thick concrete and are thermally insulated with 14-cm-thick extruded polystyrene. The thermal properties of the materials are listed in Table 2, and their optical properties are summarized in Table 3. Finally, to consider exchanges with the environment, the weather conditions were adapted to correspond to the actual geographic location of the site. The soil around the digesters is considered to be loamy, and its properties are assumed to be constant and uniform at all depths (Table 2).

Solid fresh matter and cattle slurry are mixed in a hopper to attain the 10 % total solids that are subsequently fed into the digester. The digester is fed with a mix of solid and liquid organic waste, with an annual mean composition (dry matter basis) of 45 wt% cattle slurry, 21 wt% cattle manure, 16 wt% corn stover, 6.7 wt% rye silage, 4.8 wt% Both the digester and the post-digester also independently receive a recirculated liquid digestate whose average flowrate depends on the season, with annual mean flowrates of 13.5 m 3 d -1 and 3.5 m 3 d -1 to the digester and post-digester, respectively.

Operating condition measurements and collection of meteorological data

Table 4 summarizes all of the collected operating conditions and meteorological data and their acquisition frequency. The measured quantities included the temperatures (of the digestate, biogas, and outer cover), the liquid level, and the solid organic waste and cattle slurry composition and flowrate fed into the digester and post-digester, as well as the flowrate and composition of the biogas.

A mix of poor biogas (~12 % CH 4 ) resulting from biomethane purification and raw biogas (~60 % CH 4 ) is burned in the main gas boiler (Fig. 4) for the production of heat. In 2021, the annually averaged flowrate of biomethane burned in the gas boiler was 9.5 Nm 3 CH 4 h -1 (including poor and rich biogas), with the energy efficiency of the combustion boiler being 93 %. Warm water heated in the boiler is used to heat the digester and post-digester. The annual mean temperatures of the water exiting and entering the boiler are 66 • C and 57 • C, respectively, and the thermal loss in the heat exchanger ducts is approximately 5 %.

As listed in Table 4, the required meteorological data for 2021 were collected corresponding to the geographical location of Biometharn. The daily dynamic ambient air temperature and wind speed for the site location were obtained from the Iowa Environmental Mesonet website [START_REF]Iowa Environmental Mesonet (IEM)[END_REF]. The daily values of the total precipitation (mm d -1 ) and soil temperature at the surface, together with the monthly values of the albedo, were taken from the Nasa Data Access Viewer [START_REF]Prediction of worldwide energy resources, Data acces viewer[END_REF]. These data are essential to dynamically model the temperature in the soil at depth (Section 2.2.3) and to determine the global solar irradiation on the Earth's surface under cloudy sky conditions (Sections 2.3.4 and 2.3.5).

Dynamic thermal results and experimental validation

The dynamic energy model was applied based on the specificities of the site and run for the entire period of 2021. The studied units consisted of a digester and a post-digester in series, as described in Section 3.1 and illustrated in Fig. 4, with the characteristics and operating conditions provided in Tables 123. Table 5 summarizes the modeling results, which are rearranged to represent the thermal losses and their locations per season for the entire unit (the digester and post-digester together) in 2021. The seasonal thermal consumption of the digesters ( QHE ) given in Table 5, ranging from 80 kW in summer to 140 kW in winter, correspond to the amount of thermal energy provided via the heat exchangers to maintain the temperature of the mesophilic AD in the digestate.

For ease of illustration and discussion, the dynamic modeling results of the main digester are illustrated for two contrasting months: August 2021 (Fig. 5a) and December 2021 (Fig. 5b). The operating temperature profiles in Fig. 5 are confined between the ambient temperature (T amb ) and the digestate temperature (T d ), with the exceptions of the cover, air, inner membrane, and biogas temperatures (T c , T air , T im , and T b , respectively), representing upper peaks that correspond to the peak sunshine during August (e.g., temperature values above 40 • C in Fig. 5a). This phenomenon results in cover overheating, especially during sunny periods in summer, and the energy model is capable of representing this phenomenon as well. Because of the overheating of the cover, the air, inner membrane, and biogas temperatures also reach values above 40 • C, as seen in Fig. 5a.

Both in Table 5 and Fig. 5, QHE represents the thermal energy required from the heat exchanger to maintain the digestate temperature. Peaks in the required thermal energy are observed during the night as a result of the absence of solar irradiation. The solar irradiation during the day markedly reduces the required heating power (-39 % on average), particularly in August as opposed to December.

Using Table 6, the results from the dynamic energy model and the experiments can be compared based on the dynamic measurements of the operating temperatures and the heat requirement with the measurement frequencies given in Table 4. Deviations are quantified via the root mean square error (RMSE) and normalized RMSE, with overall RSME values of approximately 1.4 • C for the temperatures and 4.2 kW (3.6 %) for the heating requirement. Fig. 6 compares the actual heat consumption of the site with the predictions made by the model on a monthly basis. The developed energy model is shown to be capable of predicting the thermal behaviors and governing temperatures, as well as the amount of thermal energy required by the heat exchanger, with acceptable errors given the scale of the unit.

Table 2

Thermal properties of the materials used in the energy model.

Table 4

List of on-site measurements and collected meteorological data and their acquisition frequency. The model-based seasonal analysis, as provided in Table 5, allows the thermal efficiency and the distribution of heat losses in the digesters and jointed gasholders up to the protective cover to be studied. The rearranged results highlight the thermal behavior of the unit depending on the season: winter (December, January, and February), spring (March, April, and May), summer (June, July, and August), or autumn (September, October, and November). The thermal losses are then graphically shown for the entire unit in Fig. 7 for the hot (summer) and cold (winter) seasons.

Regarding the overall losses of the entire unit in Table 5, for all seasons, the results indicate that the main heat losses are linked to the radiative losses from the cover to the atmosphere ( Qc,at(rad) ), the net advective digestate losses (Δ Ḣd ), and the convective cover loss to the atmosphere ( Qc,at ), representing 29 %, 22 %, and 21 %, respectively, of the total losses in winter, and 27 %, 10 %, and 26 %, respectively, of the total losses in summer (Table 5). Note that Δ Ḣd represents the net difference between the advective flows at the output and input, which decrease in summer as a result of the increased ambient temperature.

Heat losses from the cover to the atmosphere are primarily due to radiative and convective effects. The radiative losses ( Qc,at(rad) ) vary between 58 kW and 76 kW, while the solar irradiation absorbed by the cover, QSunc(rad) , varies between 42 kW in winter and 147 kW in summer. This strong variation in the solar irradiation received by the cover between the cold season (autumn and winter) and the hot season (spring and summer) is primarily due to the high level of cloudiness during the cold season. The overall radiative and convective cover losses to the atmosphere in Table 5 are smaller in winter (110 kW) than in summer (172 kW) because of the higher amount of solar radiation received in summer. The same phenomenon occurs for the overall radiative and convective wall losses to the atmosphere, with values of 18 kW and 55 kW in winter and summer, respectively.

Moreover, as shown in Table 5, the total heat loss from the liquid digestate compartment is 140 kW in winter, 123 kW in spring, 79 kW in summer, and 118 kW in autumn. Naturally, the losses are greater during winter as a result of the lower ambient temperature. Regarding the radiative and convective losses associated with the liquid digestate ( Qd,net(rad) ), winter represents the largest loss, at 37 %, because of the lower inner membrane temperature (26.3 • C in winter compared with 35.1 • C in summer) and cover temperature (11 • C in winter compared with 29.5 • C in summer), which intensifies the radiative transfer from the hot digestate to the upper gasholder zone.

Thermal losses to the underground walls and bottom ( Qd,gr ) and to the aboveground walls ( Qd,w ) represent a small component of the losses, 6 % on average, as a result of the use of efficient thermal insulation materials (e.g., extruded polystyrene) in modern anaerobic digesters and the insulating effect of the soil. The net amount of heat absorption related to the biodegradation reaction and the stirrer mixing (Δ Ḣr -ẆM ) is nearly constant over the year, remaining below 4 kW. As discussed in the previous section, losses linked to the radiative and convective heat transfer from the digestate to the cover and then to the atmosphere are among the main sources of thermal loss. Heat losses from the digestate through the membranes to the cover can be reduced via the use of multiple membranes between the digestate and the cover. We studied the addition of membranes using the developed thermal model, in which two additional equations (similar to Eqs. ( 2) and ( 4))

were included per added membrane as indicated in Section 2. The presence of insulating membranes makes it possible to reduce the convective and radiative heat received by the cover, i.e., to reduce the heat loss of the digestate. Based on an annual average, the modeling results indicate that the use of a third insulating membrane (an additional 5-cm-thick polystyrene insulation surface between the digestate and the inner flexible membrane of the gasholder) allows for a 34 % reduction in the energy requirements compared with a doublemembrane gasholder (Fig. 8). Fig. 8 illustrates that the influence of the number of insulation membranes on the total amount of heat supply required can be evaluated via energy modeling. Here, the number of membranes ranges from one (including only the protective cover) to three (including the cover, the inner membrane, and a third insulation barrier). It can be seen that switching from a single-to doublemembrane configuration provides a thermal gain of 27 %, while switching from a single-to triple-membrane system results in a total thermal gain of 52 %. b) Selective cover coatings. As shown in Table 5, the solar irradiation absorbed by the cover is significant, ranging from 42 kW in winter to 147 kW in summer as a result of the high solar absorptivity of the materials on the superior surface of the cover (α s,co,sup ). Nevertheless, a significant proportion of this energy is soon reemitted to the sky as longwave infrared radiation depending on the material emissivity (ε) and the temperature, with values ranging from 58 kW to 76 kW throughout the year.

The choice of cover materials with selective coatings may play a key role in improving the thermo-optical efficiency of the cover by reducing the longwave infrared emissivity to values lower than 0.3, while increasing the solar absorptivity to values greater than 0.95, resulting in selectivity (α s /ε) values ranging from 3 to 10, via the use of lowemissivity selective materials or paints. Based on a non-costly selective paint application with a selectivity of 3 (ε = 0.3 and α s = 0.9), the simulations indicate the following interrelated thermal loss reductions: -12 % in terms of the total heat supply ( QHE ) and -13 % in terms of the cover loss to the atmosphere ( Qc,at + Qc,at(rad) ) during winter. c) Recovery of waste heat from biogas purification. As part of the purification of biogas for injection, gas separation units such as water or amine scrubbing and Pressure Swing Adsorption (PSA) are used. In the case of high pressure compression (e.g., PSA for biomethane purification), there are compressors equipped with cooling systems allowing the heat produced by the equipment to be transmitted at a temperature above 70-80 • C to the main heating network. On an annual basis, in 2021, the heat available from PSA was estimated at 25 % of the annual mean thermal requirements for heating the digesters. d) Energy recovery from the outgoing digestate. Concentric double-tube heat exchangers, as well as spiral plate heat exchangers, can be compared with respect to the preheating of the fresh matter and slurry at the inlet using the outgoing hot digestate from the digesters following a counter-current heat exchange configuration. Heat exchanger efficiencies between 70 % and 90 % can be reasonably expected, depending on the type of heat exchanger used.

In the case of concentric double-tube heat exchangers applied to the operating conditions at Biometharn, the exchanger efficiency barely exceeds 70 % for a maximum heat exchange power of 53 kW because of the considerable length of the tubes (180 m) and the required heat exchange surface (62 m 2 ), with the internal diameters of the tubes being 10 cm and 20 cm, corresponding to a low global heat transfer coefficient (U) of 78 W m -2 • C -1 , a number of transfer units equal to 2, and ΔT ml = 8 • C.

Conversely, spiral flow heat exchangers can actually reach a 90 % or higher efficiency for the same available maximum heat exchange power with a significantly reduced heat exchange surface of 16 m 2 (an outside spiral diameter of 65 cm), corresponding to a high global heat transfer coefficient (U) of 875 W m -2 • C -1 and ΔT ml = 3 • C. Spiral flow heat exchangers represent, moreover, the advantages of compactness (i.e., a reduced device size for the required heat exchange area), a reduced rate of fouling because the curved flow path and the single flow pass limit the amount of solid particle settlement, ease of cleaning, and lower maintenance costs.

Simulations were conducted to determine the thermal gains related to the heat recovery from outgoing digestate in a biogas plant using a spiral flow heat exchanger operating at 90 % efficiency. The results indicate that the heat recovery from the exiting digestate reduces the total heat supply ( QHE ) by -44 % for digesters equipped with a doublemembrane gasholder and by -68 % for digesters equipped with a triplemembrane gasholder.

Overall, as indicated in Table 7, independent of the recovery of waste heat from the purification step, by combining the effect of a triplemembrane gasholder with that of the recovery of thermal energy from the outgoing digestate, the self-consumption of biogas is reduced to 1.4 %, on an annual basis. When the possibility of thermal recovery of waste heat from the biogas purification step exists, the anaerobic digester can achieve thermal autarky. Note that, here, the use of additional components such as spiral plate heat exchangers to recover advection heat from the outgoing digestate, and the alternative choice of construction and insulation materials as well as the alternative choices for AD design and operation have only been discussed from an energy efficiency perspective. This is because the energy and economic benefits are linked and are proportional to the energy efficiency and the reduction of the self-consumption of primary energy (biogas). For example, in the context of cogeneration, the purchase price of electricity produced from biogas benefits from an energy efficiency bonus, calculated on the basis of the electrical energy produced and the thermal energy recovered excluding the self-consumption of the plant for an energy efficiency value between 35 % and 80 %. In the context of biomethane production, the purchase price of biomethane

Table 7

Reduction of the self-consumption of biogas via multimembrane gasholders and various heat recovery strategies (annual mean). also depends on the annual production of the biogas plant. A comprehensive cost analysis is therefore a relevant subsequent component of this research that will provide additional information for decisionmaking.

Conclusions

This study presented the development of a thermal model and an operational simulation tool for semi-buried large-scale anaerobic digesters with upper multimembrane gasholders. The experimental validation, via the test case of a large-scale industrial biogas plant, demonstrated the importance of including a detailed geometric representation of the digesters and upper gasholders with multiple membrane insulation, where the longwave infrared radiative network is also fully applied to all involved surfaces. The related overall radiative and convective losses from the multimembrane gasholder represented, on average, 56 % of the overall annual thermal losses of the digester.

Seasonal analyses were used to provide an overview of the thermal performance of the digester and to make it possible to identify the main sources of digestate loss in the studied AD unit. On average, the most significant losses are the digestate advective loss (33 %), radiative digestate loss (33 %), biogas advective loss (14 %), and convective digestate loss (11 %).

Identifying the main sources of thermal loss enabled an ex ante study of the main thermal strategies for reaching thermal autarky in anaerobic digesters for large-scale biogas production. Three main strategies were proposed and investigated to evaluate the thermal autarky potential of digesters and gasholders. The effect of multimembrane insulation in the gasholder revealed that the thermal consumption of a plant could be reduced by -34 % when using a triple-membrane gasholder compared with a double-membrane gasholder. The use of a selective coating on the cover, alone, could provide a reduction of -12 % of the total heat provision, and the use of an appropriate heat exchanger recovering the advective heat of the outgoing digestate to preheat the digester feed could reduce the total heat supply by -44 % in the case of a digester equipped with a double-membrane gasholder and by -68 % in the case of a digester equipped with a triple-membrane gasholder.

Overall, the energy simulations demonstrated that the application of triple-membrane gasholders together with the use of a spiral plate heat exchanger for heat recovery from the digestate advective heat would allow the annual biogas self-consumption of a plant to be reduced to 1.4 %, i.e., only 1.4 % of the total raw biogas produced would be used for heating the plant. This result represents a cover loss reduction of -51 % and an advective digestate loss reduction of -81 % when compared with the actual operation of the site. When waste heat recovery from biogas purification is also possible and applied, the annual biogas selfconsumption of a plant can reach 0 %, i.e., achieve thermal autarky.

In addition to the useful insights drawn from the application and use of the dynamic energy model developed in this study, this model can be used as a reliable tool to conduct further studies via sensitivity analyses and optimizations of large-scale AD units, for example, studies concerning the geometry and dimensions of digesters and gasholders, scenarios of buried or semi-buried digesters, the thermal insulation of storage equipment, and the use of solar collectors for heating.

CRediT

Declaration of Competing Interest

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

  supplied by the work of the mixing systems (Thermal expansion coefficient (K -1 ) γ Split factor between the sky and the surroundings (-) δ Thickness (m) ε Emissivity (-) θ Tilt angle of the surface ( • ) μ Dynamic viscosity of the fluid (Pa s) ν Kinematic viscosity of the fluid (m 2 s -1 )

Fig. 1 .

 1 Fig. 1. General configuration of a modern wet anaerobic digester with a multimembrane gasholder.

  7) when T d > T b or following the principle of the convection on the upper side of a horizontal cold plate using the McAdams [27] correlation (Eq. (A.8)) when T d ≤ T b .

  ) and Eqs. (A.20) and (A.21), respectively, by determining the different thermal resistances involved, i.e., the wall and floor, as well as the internal and external surfaces in contact with the soil.

Fig. 3 .

 3 Fig. 3. Schematic representation of the enclosure and the view factors from an external surface to the sky, surroundings, and ground. The view factors from the walls to the sky (F w-sky ), from the walls to the surroundings at ambient air temperature (F w-sur ), and from the walls to the ground (F w-gr ) are determined according to the tilt angle of the surface θ (equal to 90 • ) and the split factor between the sky and the surroundings, γ = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 0.5(1 + cosθ) √ , based on an enclosure consisting of an exterior surface, the sky, the surroundings, and the ground surface.

Fig. 4 .

 4 Fig. 4. Schematic illustration of Site A combining a digester and a post-digester for biogas production. The dashed lines represent the water flow in the circulators heating the digester and post-digester. The red and blue lines represent the hot and cold water flows, respectively, in the main boiler ducts.

  temperature ( • C) Daily Fresh solid organic waste flowrate (kg d -1 ) Daily Cattle slurry (kg d -1 ) Daily Liquid digestate recirculated (m 3 d -1 ) Daily Power consumption of the stirrer (kWh) Daily Digestate temperature ( • C) 5 min Biogas temperature ( • C) 5 min Cover temperature ( • C) 5 min Warm water flowrate from heating system (m 3 h -1 ) 15 min Outgoing water temperature of heating system ( • C) 15 min Incoming water temperature of heating system ( • C) 15 min Biogas burned in the combustion boiler (N m 3 h -1 ) Daily Meteorological data Ambient air temperature ( • C) 30 min Wind speed (m s -1 ) 30 min Precipitation (m d -1 ) Daily Soil surface temperature ( • C) Daily Albedo Monthly Insolation clearness index KT Monthly

3. 3 .

 3 Strategies to approach thermal autarky in biogas production plants a) Multimembrane effect in the gasholder.

Fig. 5 .

 5 Fig. 5. Model predictions of the temperatures and heat supply of the digester in (a) August and (b) December 2021. The top plots show the temperature profiles: T d (light green); T b (dark green); T c (black); T w (red); T im (light blue); T air (purple); and T amb (grey). The bottom plots show the heat supply from the heat exchanger QHE .

Fig. 6 .

 6 Fig.6. Monthly model-experiment comparison for the heat supply (Q EH ) to the digesters in 2021. The heat loss from the walls to the atmosphere accounts for 9%-20% of the total losses, depending on the season. As discussed above, this relatively smaller contribution, when compared with the cover, is primarily explained by the thickness of the walls, as well as the use of thick extruded polystyrene as insulation over the external surface of the digester walls.

Fig. 7 .

 7 Fig. 7. Distribution of the overall heat losses for the entire biogas plant (i.e., the digester, post-digester, and gasholders) in the winter and summer of 2021.

Fig. 8 .

 8 Fig. 8. Reduction in the total heat supply for the multimembrane gasholder configurations (annual mean, 2021).

  authorship contribution statement M. Avila-Lopez: Writingoriginal draft, Methodology, Software, Investigation. C. Robles-Rodriguez: Resources. L. Tiruta-Barna: Resources, Writingreview & editing. A. Ahmadi: Conceptualization, Methodology, Software, Investigation, Writingreview & editing.

Table 1

 1 Main characteristics of the digester and post-digester. , and greases, 3 wt% bract, and 1.2 % other culture waste. The digester is fed with an annual mean flowrate of 24.6 t d -1 of fresh solids and 20.1 t d -1 of cattle slurry. The post-digester is directly fed by the outgoing digestate and biogas coming from the digester.

	Characteristics	Digester	Post-
			digester
	Digester diameter (m)	25	18
	Digester height (m)	6
	Depth of the buried part of the digester (m) Volume (m 3 ) Confined biogas volume (m 3 )	2.5 2950 310	1530 160
	Tilt angle of the cover ( • )	22.8	31.3
	Concrete wall thickness and bottom (m)	0.2
	Thickness of the insulation membrane -ethylene	0.005
	propylene diene monomer (EPDM) (m)
	Thickness of the cover membrane -polyvinyl chloride	0.005
	(PVC) (m)
	Exterior wall insulation thickness -extruded polystyrene	0.14
	(XPS) (m)
	Number of membranes	2
	Height between the digestate and the cover (m)	2.52
	Mechanical stirring system (kW) Annual mean fresh matter flowrate to digester (t d -1 ) Annual mean slurry flowrate to digester (kg d -1 )	15 24.6 20.1	4
	Total solids in fresh matter (%)	30-35
	Total solids in feed to digester (%) Annual mean raw biogas flowrate (Nm 3 h -1 )	10 295.5
	Annual mean raw biogas CH 4 composition (%vol)	58.5
	Annual mean raw biogas CO 2 composition (%vol)	41.3
	Type of soil	Loam

Table 3

 3 Optical properties of the construction and insulation materials.

	Thermal	Concrete	XPS	EPDM	PVC	Loam soil
	property
	Thermal	1.4	0.034	0.25	0.17	2.2
	conductivity k(W m -1 K -1 )
	Heat capacity cp(J kg -1 K -1 )	880	1450	1000	900	1315.8
	Density ρ(kg m -3 )	2300	42.5	1150	1390	1710
	Thermal diffusivity α(m 2 s -1 )	6.92 × 10 -7	5.52 × 10 -7	2.17 × 10 -7	1.36 × 10 -7	9.78 × 10 -7
	Thermal property	Concrete	XPS	EPDM	PVC
	Absorptivity α	0.94		0.86	0.86	0.86
	Emissivity ε	0.94		0.86	0.86	0.86
	Solar absorptivity αs	0.60		0.74	0.74	0.74

Table 5

 5 Seasonal distribution of the absolute digestate and the overall heat flows and thermal losses for the entire biogas production plant (i.e., digester, post-digester, and gasholders) in 2021.

			Winter	Spring	Summer	Autumn
	Flowrates and operating temperatures
	Incoming fresh	ṁf (td -1 )	55.0	59.5	59.1	59.2
	matter
	Produced biogas	ṁb (Nm 3 h -1 )	301.9	299.4	294.6	285.7
	Digester	T d ( • C)	40.1	40.9	40.9	40.8
	temperature
	Post-digester	T post-d ( • C)	40.8	41.0	40.9	40.8
	temperature
	Overall energy input
	Heat supply	QHE (kW)	139.7	122.5	79.3	117.7
	required from the
	heat exchanger
	Heat dissipated by	ẆM(kW)	17.1	17.1	17.1	17.1
	the stirrer
	Shortwave solar irradiance	QSun-c(rad) (kW)	41.7	116.1	147.0	68.9
	absorbed by the
	cover
	Shortwave solar irradiance	QSun-w(rad) (kW)	14.9	42.1	53.4	24.6
	absorbed by the
	walls
	Digestate energy loss
	Digestate advective	Δ Ḣd (kW)	42.4	42.6	26.7	36.9
	loss
	Biogas leaving the	Δ Ḣd,b (kW)	14.6	15.8	17.0	14.9
	digestate
	advective loss
	Convective heat	Qd,b(rad) (kW)	19.2	13.1	5.7	13.7
	transfer from
	digestate to
	biogas
	Convective/	Qd,w (kW)	2.7	2.0	1.0	2.0
	conductive heat
	transfer from
	digestate to
	aboveground
	walls
	Conductive heat	Qd,gr (kW)	5.5	4.9	3.4	4.5
	transfer from
	digestate to
	underground
	walls and ground
	Net radiative flow	Qd,net(rad) (kW)	52.3	40.1	22.1	41.1
	from the digestate
	surface
	Net heat of	Δ Ḣr -ẆM(kW)	3.1	4.2	3.6	4.5
	biodegradation
	and mixing
	Total loss	(kW)	139.8	122.7	79.5	117.6
	Overall energy loss
	Digestate advective	Δ Ḣd (kW)	42.4	42.6	26.7	36.9
	loss
	Biogas advective	Δ Ḣb (kW)	13.7	15.0	16.6	14.1
	loss
	Air ventilation	Δ Ḣair(kW)	3.3	3.7	3.2	3.0
	advective loss
	Conductive heat	Qd,gr (kW)	5.5	4.9	3.4	4.5
	transfer from
	digestate to
	underground
	walls and ground
	Convective heat	Qw,at (kW)	6.2	22.4	24.9	10.8
	transfer from the
	walls to the
	ambient air
	Convective transfer	Qc,at (kW)	41.7	78.5	71.4	42.5
	from the cover to
	the ambient air
	Convective heat	Qc,rain (kW)	10.3	14.8	24.0	17.7
	loss due to rain
		Qw,at(rad) (kW)	12.3	22.4	29.7	16.5

Table 5

 5 (continued )

	Winter	Spring	Summer	Autumn

Table 6

 6 Monthly averages of the main results from the dynamic modeling in 2021 compared with on-site experimental measurements in terms of the heat supply and operating temperatures (the root mean square error (RMSE) and normalized RMSE results are shown in gray).

		Total	Digester				Post-digester
		Q HE (kW)	Td( • C)	Tb( • C)	Tw( • C)	Td( • C)	Td( • C)	Tb( • C)	Tw( • C)	Td( • C)
	January	153.1	40.0	32.8	5.5	8.6	40.9	33.0	5.6	8.2
	February	131.2	40.1	34.1	12.0	13.9	41.2	34.4	12.0	13.7
	March	132.6	40.9	35.1	13.5	16.0	41.3	34.9	13.4	15.5
	April	129.4	41.0	35.8	16.8	19.0	40.9	35.3	16.8	18.3
	May	105.5	41.0	36.9	22.6	23.8	40.9	36.4	22.5	23.1
	June	78.9	41.1	38.2	28.9	29.3	41.0	37.8	28.8	28.6
	July	77.8	41.0	38.5	30.2	30.6	41.0	38.1	30.0	29.5
	August	81.2	40.7	38.0	29.0	29.6	40.7	37.7	28.9	29.0
	September	94.8	40.8	37.1	24.8	25.2	40.8	36.8	24.8	24.8
	October	116.7	40.9	35.7	16.6	18.6	40.9	35.3	16.6	18.1
	November	141.5	40.7	34.1	9.3	12.0	40.9	33.8	9.3	11.5
	December	134.7	40.3	33.5	8.3	10.9	40.3	33.0	8.3	10.5
	Year average	114.8	40.7	35.8	18.1	19.8	40.9	35.6	18.1	19.2
	RMSE	4.2 kW	1.38 • C
	Normalized RMSE	3.63 %	6.2 % (

• C basis)

Acknowledgements

The authors acknowledge GRDF (Gaz Réseau Distribution France) within the framework of the Biogas Innovation Chair established at INSA Toulouse, France, for the financial support given to the ECOFEV project. The authors also thank BIOMETHARN (an industrial biogas plant in France) for agreeing to provide real-time data of their operation for one year. We thank Martha Evonuk, PhD, from Evonuk Scientific Editing (http://evonukscientificediting.com) for editing a draft of this manuscript.

Data availability

The data that has been used is confidential.

Appendix

Appendix A. . Convective and conductive heat transfer formulations A.1 Dimensionless numbers in heat transfer phenomena

Non-dimensional numbers, primarily the Nusselt number (Nu), Prandtl number (Pr), Grashof number (Gr), Rayleigh number (Ra), and Reynolds number (Re), are used in the thermal modeling of the different conductive and convective heat transfer terms involved in the balance equations. These numbers are briefly described below.

The Nusselt number is the ratio of convection to pure conduction heat transfer at a boundary in a fluid. Determining the Nusselt number is a key step in estimating the heat transfer coefficient h.

Most correlations in the literature concerning the Nusselt number for convection are expressed in terms of Nu = f(Ra, Pr) or Nu = f(Re, Pr) for natural or forced convection, respectively.

The Prandtl number is described as the ratio of the momentum to thermal diffusivities.

The Grashof number measures the ratio of the buoyancy forces to viscous forces; this number plays a key role in natural convection.

For ideal gases, β, or the thermal expansion coefficient, is taken to be 1/T i . The Rayleigh number is the product of the Grashof and Prandtl numbers.

The Reynolds number represents the ratio of the inertial to viscous forces and controls the flow patterns of fluids (i.e., laminar, transitional, or turbulent); it is primarily used in forced convection correlations.

A.2 Natural convection correlations

A.2.1 Horizontal plates

For horizontal plates, the flow patterns and heat transfer depend strongly on whether the surface is cold or hot and on whether it is facing upward or downward. For the upper surface of a hot plate or the lower surface of a cold plate, the Nusselt number is calculated using the following equations [START_REF] Mcadams | Heat Transmission[END_REF].

Pr ≥ 0.7 (A.6)

All Pr (A.7).

For the lower surface of a hot plate or the upper surface of a cold plate, the Nusselt number is determined via Eq. (A.8) [START_REF] Mcadams | Heat Transmission[END_REF].

Pr ≥ 0.7 (A.8)

In both cases, the characteristic length used in the Rayleigh number for a circular horizontal plate is half the plate diameter.

A.2.2 Vertical plates

The Nusselt number for natural convection in vertical plates can be obtained using the following correlations [START_REF] Churchill | Correlating equations for laminar and turbulent free convection from a vertical plate[END_REF]: (A.10)

Here, the characteristic length used in the Rayleigh number is the length of the vertical plate.

A.3 Forced convection correlations

A.3.1 Flat plates in parallel flow

The Nusselt number for the forced convection of a parallel flow over a flat plate depends on the Reynolds number and is calculated using the following correlations [START_REF] Incropera | Fundamentals of heat and mass transfer[END_REF].

A.3.2 Cylinder in a cross flow

The Nusselt number for the forced convection of fluid motion normal to the axis of a circular cylinder can be obtained using the correlation proposed by Churchill and Bernstein [START_REF] Churchill | A correlation equation for forced convection from gases and liquids to a circular cylinder in crossflow[END_REF], which is an equation that covers the entire Re and Pr ranges. This correlation has the following form:

For (Re × Pr) ≥ 0.2 (A.13).

Here, the corresponding characteristic length is the diameter of the cylinder.

A.4 Thermal losses to the ground using the EN ISO 13,370 standard

The global coefficient of the heat transfer to the ground is determined by separating the losses through the underground wall to the upper soil layer and through the bottom floor to the soil at a depth z.

For the digester floor,

where, w is the thickness of the adjacent walls, R si the internal surface thermal resistance, R floor the floor thermal resistance, R se the external surface thermal resistance, k soil , the thermal conductivity of the soil With, R w the wall thermal resistance.

Appendix B. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.fuel.2022.126978.