Ollivier Laccourreye
email: ollivier.laccourreye@aphp.fr

Haitham Mirghani

Diane Evrard

Pierre Bonnefont

Lydia Brugel

Frederic Tankere

André Coste

Beatrix Barry

Bertrand Baujat

Sarah Atallah

Romain Kania

Impact of the first month of Covid 19 lockdown on oncologic surgical activity in the Ile de France region university hospital otorhinolaryngology departments

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Abstract. Objective: To evaluate the impact of the first month of lockdown related to the Covid 19 epidemic on the oncologic surgical activity in the Ile de France region university hospital otorhinolaryngology departments. Material and Methods: A multicenter prospective observational assessment was conducted in 6 university hospital otorhinolaryngology departments (Paris Centre, Nord, Est and Sorbonne) during the 1-month periods before (Month A) and after (Month B) lockdown on March 17, 2020. The main goal was to evaluate lockdown impact on oncologic surgical activity in the departments. Secondary goals were to report population characteristics, surgery conditions, postoperative course, progression of Covid status in patients and surgeons, and adverse events. Results: 224 procedures were performed. There was 10.9 % reduction in overall activity, without significant difference between departments. Squamous cell carcinoma and larynx, hypopharynx, oropharynx, oral cavity and nasal cavity and sinus locations were predominant, at 79% and 75.8 % of cases respectively, with no significant differences between months. T3/4 and N2/3 tumors were more frequent in Month B (p = .002 and .0004). There was no significant difference between months regarding surgical approach, type of reconstruction, postoperative course, tracheotomy and

Introduction.

In France, the first three cases of the Covid-19 epidemic, in patients of Chinese origin who had been staying in Wuhan, were diagnosed on January 24, 2020 (https://fr.m.wikipedia.org).

On March 17, lockdown was implemented to combat the epidemic, and on March 25 the French Society of ENT (SFORL) published guidelines for head and neck cancer surgery [START_REF] Fakhry | French consensus on management of head and neck cancer surgery during COVID-19 pandemic[END_REF].

The main aim of the present study was to assess the impact of this month of lockdown and of the SFORL guidelines in university hospital otorhinolaryngology (ORL) departments in the Ile de France region, where the epidemic was at its most severe. Secondary objectives were to report population characteristics, surgery conditions, postoperative course, progression of Covid status in patients and surgeons, and adverse events.

Material and Method.

Six university hospital ORL departments in the Ile de France region (Paris Centre, Paris Nord, Paris Est and Paris Sorbonne) practicing oncologic surgery took part in a multicenter prospective observational study.

The main aim of the present study was to assess the impact of this month of lockdown and of the SFORL guidelines in university hospital otorhinolaryngology (ORL) departments in the Ile de France region, where the epidemic was at its most severe. Secondary objectives were to report population characteristics, surgery conditions, postoperative course, progression of Covid status in patients and surgeons, and adverse events. Review board approval was not sought, as the study did not come under French legislation for studies

Results.

1.Activity.

224 oncologic procedures were performed, with a 10.9% decrease between Month A (118) and Month B (106). Progression varied between departments, with 15.7%, 24.2% and 23.3% decrease in 3 cases and 4.3%, 20% and 50% increase in the other 3 (Appendix 1), without significant difference (p= 0.76).

Appendix 1 shows diagnostic procedures, with a 6.3% decrease (Month A: 47; Month B: 44). Progression varied between departments, with 14.2%, 33.3% and 50% decrease in 3 cases and 20% and 22.2% in 2; the other department had no diagnostic activity during the study period (Appendix 1). Again these differences were non-significant (p = 0.81).

Appendix 1 shows therapeutic procedures (n=133), with a 12.6% decrease (Month A: 71 -Month B: 62). Progression varied between departments, with 7.6%, 19%, 31.2% and 50% decrease in 4 cases and 20% and 83.3% increase in 2 (Appendix 1), without significant difference (p = 0.54).

Population, approaches, procedures and postoperative course.

Table 1 shows the main population characteristics. T3/4 and N2/3 tumor was significantly more frequent in Month B (Table 1). Squamous cell carcinoma and larynx, hypopharynx, oropharynx, oral cavity, nasal cavity and facial sinus locations predominated, at respectively 79% (177/224) and 75.8% (170/224) of procedures (Appendix 1), without significant inter-month difference (p respectively 0.74 and 0.53).

There were no significant inter-month differences in surgical approach, lymph-node surgery, type of reconstruction, tracheotomy and permanent tracheostomy or in postoperative course in terms of complications and mortality (Table 2 -Appendix 1). Type of mask differed significantly (p < 0.0001), with an increase in use of FFP2 masks from 11% (13/118) in Month A to 66.2% (66/106) in Month B. There were no significant differences in tracheotomy time or nasogastric feeding, intensive care unit and department stay or hospital stay (Table 3, Appendix 1).

Covid status, adverse events.

Percentage nasopharyngeal swab and polymerase chain reaction (PCR) analysis on the eve of surgery varied significantly between months and departments (p < 0.0001): from 10.1% in Month A to 33% in Month B (Table 1) and 5.6%, 10%, 10.3%, 17.1%, 18.1% and 40,4% according to department (Appendix 1). No Month A patients were Covid-19-positive, versus 3 in Month B (Appendix 1). The first of these 3 was asymptomatic, scheduled for petrectomy and parotidectomy, and had negative PCR on eve of surgery. At day 6, the surgeon was informed that the examination was in fact slightly positive; PCR at discharge was negative. The patient was readmitted on day 17 for local superinfection, with PCR again positive and becoming negative during the hospital stay. The second patient had suggestive Covid-19 symptomatology and positive PCR at eve of surgery scheduled for glossectomy with free flap and neck dissection, which was therefore postponed for 2 weeks. PCR was still slightly positive at eve of surgery, but the operation was maintained due to increased tumor volume and absence of respiratory symptoms. The third patient had moderate pulmonary pathology and underwent salvage total pharyngolaryngectomy with free flap after a 2-week postponement. The other 291 patients did not develop signs of Covid-19 during the study period or following month; nor did the senior and junior physicians involved in surgery.

Discussion.

The prospective study in more than 500 Chinese hospitals published in Lancet Oncology in March 2020 was the first analysis of the relation between Covid-19 and cancer, and suggested that persons with or treated for cancer (mainly lung or colorectal) showed extra risk of viral infection and of severe disease, especially in case of chemotherapy or surgery in the month preceding infection [START_REF] Liang | Cancer patients in SARS-CoV-2 infection: a nationwide analysis in China[END_REF]. In China, these findings led to guidelines advocating changing certain practices during the epidemic in digestive and lung cancer patients, with 2 main aims in view: to limit high-risk situations such as surgery, and to limit contact, especially with health-care establishments. At the same time, several societies, groups and associations specializing in head and neck cancer surgery put forward various advice [START_REF] Fakhry | French consensus on management of head and neck cancer surgery during COVID-19 pandemic[END_REF][START_REF] Cai | Treating head and neck tumors during the SARS-CoV-2 epidemic, 2019-2020: Sichuan cancer hospital[END_REF][START_REF] Ansarin | Surgical management of head and neck tumours during the SARS-CoV (Covid-19) pandemic[END_REF][START_REF] Saibene | Role and management of a head and neck department during the COVID-19 outbreak in Lombardy[END_REF][START_REF] Day | Head and neck oncology during the COVID-19 pandemic: Reconsidering traditional treatment paradigms in light of new surgical and other multilevel risks[END_REF][START_REF] Givi | Safety recommendations for evaluation and surgery of the head and neck during the COVID-19 pandemic[END_REF][START_REF] Kulcsar | High risk of COVID-19 infection for head and neck surgeons[END_REF][START_REF] Rassekh | Management of the difficult airway in the COVID-19 pandemic: Illustrative complex head and neck cancer scenario[END_REF][START_REF] Shanti | Considerations in the evaluation and management of oral potentially malignant disorders during the COVID-19 pandemic[END_REF][START_REF] Yuen | Head and neck cancer care in the COVID-19 pandemic: A brief update[END_REF][START_REF] Kowalski | COVID-19 pandemic: Effects and evidence-based recommendations for otolaryngology and head and neck surgery practice[END_REF][START_REF]Position statement: otolaryngologists and the COVID-19 pandemic[END_REF]. The first idea was a form of triage, adapting the interval to surgery according to the potential danger of the cancer. The second was to reduce head and neck oncologic surgery activity so as to be able to concentrate anesthesiology personnel and equipment on intensive care units that were experiencing exponential rises in admissions due to Covid-19. The third, given a threatened shortage of protective equipment (gloves, masks, etc.), was to protect ENT physicians against possible infection, ENT and head and neck surgery being at high risk of transmission due to tracheotomy and mucosal resection near infected areas.

By comparing data for the months preceding and following the start of lockdown on March 17 from six hospital departments, the present study reflects the impact of lockdown and of the SFORL guidelines [START_REF] Fakhry | French consensus on management of head and neck cancer surgery during COVID-19 pandemic[END_REF] on the practice of head and neck oncologic surgery in university hospital departments in the Ile de France region, where the Covid-19 epidemic was most rife. The first finding was that head and neck oncologic surgery was relatively stable in these departments over the two months straddling the start of lockdown, with an overall fall of 10.9%. The decrease was 6.3% for diagnostic procedures and 12.6% for therapeutic procedures, without significant differences between the two time periods or the six departments; nor did it involve any particular change in clinical profile, with 79% squamous cell carcinoma (Table 1) located for 75.8% in the larynx, hypopharynx, oropharynx, oral cavity or nasal cavities and facial sinuses (Appendix 1). These data, and the fact that T3/4 and N2/3 tumors were significantly more frequent in Month B (Table 1) testify in our opinion to the effort made by surgery teams to avoid any adverse impact on survival of delaying treatment for these tumors [19][20][START_REF] Schutte | Impact of time to diagnosis and treatment in head and neck cancer: a systematic review[END_REF], given that most head and neck cancers double in volume within 1 to 3 months [START_REF] Jensen | Tumor progression in waiting time for radiotherapy in head and neck cancer[END_REF] and that the first month of lockdown in France witnessed a dramatic fall in community ENT consultations, which are where the vast majority of such tumors are first detected; community ENT physician in Réunion Island reported a mean decrease of 75.6% in their consultations over this period [START_REF] Rubin | Impact de l'épidémie SARS-CoV-2 sur les consultations libérales d'otorhinolaryngologie lors du premier mois de confinement de l'île de la Réunion en 2020[END_REF].

Months A and B did not significantly differ in surgical approach, neck dissection, type of reconstruction, postoperative course, complications rate and mortality, intensive care stay, tracheotomy time and nasogastric feeding time, or hospital stay (Table 3 -Appendix 1), confirming the commitment of the surgery teams participating in this prospective study, who did not fundamentally alter their management of cancer and patients during the epidemic. Secondly, although two studies [START_REF] Rampinelli | Reshaping head and neck reconstruction policy during the COVID-19 pandemic peak: experience in a frontline institution[END_REF][START_REF] Patel | Early institutional head and neck oncologic and microvascular surgery practice patterns across the United States during the SARS-CoV-2 (COVID19) pandemic[END_REF] reported a decrease in the use of micro-anastomosed free flaps for reconstruction, the present data point to no adverse effects (Table 2). And a third finding of interest was the absence of symptoms suggestive of Covid-19 infection over the study period and following month in the 221 patients uninfected at the beginning of the study period and in the senior and junior surgeons involved in the 224 procedures, testifying in our opinion to good management of the epidemic within the hospitals of the participating departments. It is notable in this regard that a lack of protective masks, feared at one point, did not in fact affect the head and neck surgery departments and that the significant change (p < 0.0001) in the type of mask used by the surgeons in the 224 procedures, with FFP2 mask use increasing from 11% in Month A to 66.2% in Month B, was due not to equipment shortage but to individual choice, as seen in the significant variation (p < 0.0001) in type of mask when comparing the six departments.

These positive findings should be weighed against the fact that only 1.3% of the patients operated on during the study period (n=3) were affected by the virus. The treatment these 3 patients received, detailed in the result section, was in line with the SFORL guidelines [START_REF] Fakhry | French consensus on management of head and neck cancer surgery during COVID-19 pandemic[END_REF], which recommend deferring surgery in Covid-19-positive patients and adapting the delay to the potential gravity of the cancer, and allowed for successful management of the tumors without any adverse event.

Conclusion.

The present data highlight that the Covid-19 epidemic and March 17 lockdown only very slightly impacted the cancer surgery care pathway in the university hospital ORL departments in the Ile de France region, protecting patients' chances of survival, without spreading the virus. The data support the decentralized strategy operating in the hospitals involved in the study, and validate the SFORL guidelines. They suggest that surgical management of head and neck cancer should be pursued along the same lines in case of a second wave of the epidemic in Ile de France.

  nasogastric feeding-tube time, intensive care stay or hospital stay. None of the Month A patients were Covid-19-positive, versus 3 in Month B, without adverse events. None of the otorhinolaryngologists involved in the procedures developed symptoms suggesting Covid-19 infestation. Conclusion: The present study underscored the limited impact of the Covid-19 epidemic and lockdown on surgical diagnosis and cancer surgery in the Ile de France university otorhinolaryngology departments, maintaining chances for optimal survival without spreading the virus. Key-words: Covid-19, cancer, otorhinolaryngology.

  involving human subjects [2]. An Excel database (Appendix 1) available to each department to collate data on procedures performed during the month preceding (Month A: February 17 to March 17, 2020) and the month following (Month B: March 18 to April 18) initiation of lockdown. Procedures comprised: procedures resulting in diagnosis of cancer (endoscopy-biopsy, tonsillectomy, adenectomy) and cancer resection and treatment of acute complications (tracheotomy for tumoral dyspnea, salivary leakage, radionecrosis). Diagnostic procedures (endoscopy) not resulting in diagnosis of cancer, functional surgery performed after cancer treatment (tracheostomy enlargement, dilation of post-treatment stenosis, implantation or change of vocal prosthesis, tracheotomy and laser treatment for post-treatment dyspnea, reconstructive surgery) as well as head and neck surgeries for benign tumor, precancerous pathology (dysplasia, in-situ carcinoma) or basal cell carcinoma were excluded. Data for included oncologic procedures comprised: patient data (age, gender, body-mass index, smoking status, Charlson comorbidity index [3]), tumor data (location, histology, T, N and M stage [4]), diagnostic or therapeutic, surgical approach (transcutaneous, transmucosal or mixed), type of lymph-node surgery, laryngotracheal axis treatment (laryngeal disobstruction, tracheotomy, permanent tracheostomy), postoperative course (complications, intensive care unit and department stay, tracheotomy time, nasogastric feeding tube duration, hospital stay), and type of surgical mask used. In Month B, data further comprised adverse events (modification or implementation of emergency surgery, unscheduled admission, major complications and/or Covid-19-related death) and Covid-19 status progression in patients and surgeons during the study period and the following month. Inter-month comparison used Chi2 or Fisher test for qualitative variables and nonparametric Mann-Whitney U test for quantitative variables, on StatView software (StatView, USA). The significance threshold was set at p= 0.005, in line with improved scientific method [5,6].

Table 1 :

 1 Main clinical characteristics of the study population (M: median; RT: radiation therapy; H&N: head and neck; SCC: squamous cell carcinoma; *: distribution based on therapeutic procedures).

		Global	Month A	Month B	P
	Gender (female/male)	65 / 159	38 / 80	27 / 79	0.30
	Age (years)	18-92 -M:63	18-92 -M:64	29-85 -M:63	0.92
	Body-mass index	13-39 -M:24	13-37 -M:23	13-39 -M:25	0.037
	Smoking
	Pack-years	0/150 -M:28	0/100 -M:26	0/150 -M:30	0.46
	Active (yes/no)	68/50	34 / 80	34 / 70	0.66
	Comorbidity
	Charlson index [2]	0-10 -M:2	0-10 -M:2	0-8 -M:2	0.26
	H&N cancer (yes/no)	46/178	31/87	15/ 91	0.031
	Non-H&N cancer (yes/no)	39/185	20/96	19/87	0.86
	Cervical RT (yes/no)	18/206	11/107	7/99	0.62
	Type of cancer
	SCC (yes/no)	177/47	99/26	85/21	0.74
	T (1 and 2/3 and 4)*	76/38	48/12	28/26	0.002
	N (0 and 1/2 and 3)*	103/27	63/6	40/21	0.0004
	M (0/1)*	124/6	68/1	56/5	0.098
	Initial Covid status
	Not tested/negative/positive	177/44/3	106/12/ 0	71/32/ 3	< 0.001

Table 2 :

 2 Procedures (PB: parotid biopsy; UND: unilateral neck dissection; TM: transmucosal; TC: transcutaneous; MX: mixed transmucosal-transcutaneous; BND: bilateral neck dissection; TRS: permanent tracheostomy; TRT: tracheotomy; LD: laryngeal disobstruction; MAF: micro-anastomosed flap, PMF: pediculated muscle flap; LF: local flap; FF: free flap).

		Global	Month A	Month B	P
	Diagnostic
	Endoscopy/tonsillectomy	81/3	44/1	37/2	0.59
	Adenectomy/UND/BP	7/2/1	2/1/1	5/1/0	-
	Therapeutic
	Tumoral
	Approach (TO/TM/MX)	43/29/39	24/16/21	19/13/18	0.033
	Trachea (TRS/TRT/LD)	28/18/1	14/9/0	14/9/1	0.61
	Nodal
	Yes/non	85/48	43/28	42/20	0.46
	Adenectomy/UND/BND	14/26/45	7/11/25	7/15/20	0.56
	Reconstruction
	Yes/no	35/98	16/55	19/43	0.32
	LF/PMF/MAF/FF	12/10/9/3	4/5/5/2	8/5/4/1	0.64

Table 3 :

 3 Postoperative course (M: median).

		Global	Month A Month B	P
	Postoperative course
	Simple/complicated	33/191	17/101	16/90	0.99
	ENT complications (yes/no)	21/203	11/107	10/96	0.99
	General complications (yes/no)	16/208	8/110	8/98	0.99
	Death (yes/no)	3/221	2/116	1/105	0.99
	Times (days)
	Intensive care service

& neck oncology! Eur Ann Otorhinolaryngol Head Neck Dis 2019 Dec 13. pii: S1879-7296(19)30210-8. doi: 10.1016/j.anorl.2019.12.007

Acknowledgments. The authors thank the Progrès 2000 Association and Mr Clément

Jourdaine for technical support, and Ms Sophie Guiquerro for precious help with the literature search.

Disclosure of interest:

The authors have no conflicts of interest to disclose.


Robert Caulk
email: rob.caulk@gmail.com

Luc Sholt Ès

Marek Krzaczek

Bruno Chareyre

A pore-scale Thermo-Hydro-Mechanical coupled model for particulate systems

Keywords: DEM, pore network, heat transfer, fluid flow, thermo-hydro-mechanical couplings 1 Introduction

A pore scale numerical method dedicated to the simulation of heat transfer and associated thermohydro-mechanical couplings in granular media is described. The proposed thermo-hydro-mechanical approach is based on an existing hydro-mechanical model that combines the discrete element method for simulating the mechanical behavior of dense sphere packings with the finite volume method for simulating pore space fluid flow and the hydro-mechanical coupling. Within the hydro-mechanical framework, the pore space is discretized as a tetrahedral network defined by the triangulation of discrete element method (DEM) particle centers. It is this discretization of DEM particle contacts and tetrahedral pore spaces that enables the efficient conductive and advective heat transfer models proposed herein. In particular, conductive heat transfer is modeled explicitly between and within solid and fluid phases: across DEM particle contacts, between adjacent tetrahedral pores, and between pores and incident particles. Meanwhile, advective heat transfer is added to the existing implicit fluid flow scheme by estimating mass energy flux from pressure induced fluid fluxes. In addition to the heat transfer model, a thermo-mechanical coupling is implemented by considering volume changes based on the thermal expansion of particles and fluid. The conduction and advection models are verified by presenting comparisons to an analytical solution for conduction and a fully resolved numerical solution for conduction and advection. Finally, the relevance of the fully coupled thermo-hydro-mechanical model is illustrated by simulating an experiment where a saturated porous rock sample is subjected to a cyclic temperature loading.

schemes for THM computations rely on finite difference approximations of the equations governing fluid flow and heat transfers (see, e.g., [START_REF] Al-Arkawazi | Modeling the heat transfer between fluid-granular medium[END_REF]; [START_REF] Kloss | Models, algorithms and validation for opensource DEM and CFD-DEM[END_REF]; [START_REF] Shimizu | Three-dimensional simulation using fixed coarse-grid thermal-fluid scheme and conduction heat transfer scheme in distinct element method[END_REF]; Zhou et al. (2009)). The interaction between the solid and fluid phases rely on empirical closures provided in the form of correlations required to depict the momentum exchange as well as heat and mass transfers [START_REF] Deen | Review of discrete particle modeling of fluidized beds[END_REF]). These closure correlations are often empirical but can also be derived from direct numerical simulations (DNS) [START_REF] Kruggel-Emden | Direct numerical simulation of coupled fluid flow and heat transfer for single particles and particle packings by a lbm-approach[END_REF]. Although DNS-DEM models can be used to study THM processes in particulate systems [START_REF] Deen | Direct numerical simulation of flow and heat transfer in dense fluid-particle systems[END_REF]), the small mesh size to particle size ratio results in high computational effort. Thus, DNS-DEM models are restricted to systems comprised of a smaller number of particles than CFD. An alternative particulate THM coupling is based on the lattice Boltzmann method (LBM) [START_REF] Yang | Direct numerical simulations of particle sedimentation with heat transfer using the lattice boltzmann method[END_REF]; [START_REF] Zhang | Particulate immersed boundary method for complex fluid-particle interaction problems with heat transfer[END_REF]). However, the coupling strategies rely on computations of distribution functions that require an accurate representation of solid-fluid boundaries which can be both numerically and computationally challenging [START_REF] Peng | Acomparative study of immersed-boundary and interpolated bounce-back methods in lbe[END_REF]).

An example of a DEM based THM model, presented by [START_REF] Tomac | Formulation and implementation of coupled forced heat convection and heat conduction in DEM[END_REF], includes convective and conductive heat transfer processes in 2-D. Although the scheme neglects heat induced fluid expansion and heat conduction within the fluid phase, it applies a well tested pipe network model initially proposed by [START_REF] Cheng | Evaluation of effective thermal conductivity from the structure of a packed bed[END_REF] and thoroughly verified by [START_REF] Wanne | Bonded-particle modeling of thermally fractured granite[END_REF] and [START_REF] Feng | Discrete thermal element modelling of heat conduction in particle systems: Pipe-network model and transient analysis[END_REF].

The 3-D THM coupled model presented here is based on the framework of the pore-scale finite volume (PFV) scheme initially proposed by [START_REF] Chareyre | Pore-Scale Modeling of Viscous Flow and Induced Forces in Dense Sphere Packings[END_REF] for up-scaling incompressible viscous flow and later extended to compressible flow by [START_REF] Scholtès | Modeling wave-induced pore pressure and effective stress in a granular seabed[END_REF]. The scheme is derived from the pore network (PN) models, which are applied widely for simulating a variety of porous media processes [START_REF] Blunt | Flow in porous media -pore-network models and multiphase flow[END_REF]. The model is implemented in the Yade DEM open source software ( Šmilauer V. et al., 2015) and is oriented toward dense grain packing applications as encountered in geomechanics. The proposed THM scheme combines four heat transfer models: a particle-particle conduction model, a particle-fluid conduction model, a fluid-fluid conduction model and a heat advection model. In addition, thermo-mechanical couplings are considered through both effects of fluid thermal expansion and particle thermal expansion. In summary, a set of equations governing the hydraulic and thermal schemes are presented for the derivation of the geometrical considerations and numerical couplings. Next, a validation exercise is provided where each component of the proposed THM model is challenged. First, the solid-solid conduction scheme is verified against the 1-D analytical solution of the classic heat conduction equation.

Second, the convective heat transfer model (encompassing conduction combined with advection) is compared to a fully resolved CFD solution considering the flow of a hot fluid through a cold particle assembly. Third, the full THM scheme is used to simulate an experiment where a saturated rock sample is subjected to thermal loading.

Methods

Mechanical scheme

The Lagrangian discrete element method (DEM) represents the mechanical behavior of a particulate system as a collection of interacting masses, where interactions between masses follow predefined forcedisplacement laws. Similar to the original DEM scheme [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], the present scheme implemented in Yade open DEM integrates particle positions through time according to Newton's second law of motion, which can be written for the whole system as:

Mẍ = f (1)
with ẍ the vector containing each particle acceleration, M the diagonal matrix of particle masses, and f the vector containing the total forces applied on the particles. The explicit central finite difference time stepping scheme integrates the particle acceleration from the current step to update the particle position at the next step (see ( Šmilauer V. et al., 2015) for details of the implementation). The inter-particle forces, f i j , depend on a contact model, F i j , such that:

∂ f i j ∂t = F i j (x i , x j , ẋi , ẋ j ) (2)

Compressible flow scheme

The pore-finite volume (PFV) scheme is used to model compressible fluid flow between solid particles [START_REF] Scholtès | Modeling wave-induced pore pressure and effective stress in a granular seabed[END_REF]. As presented in [START_REF] Chareyre | Pore-Scale Modeling of Viscous Flow and Induced Forces in Dense Sphere Packings[END_REF], the PFV-DEM coupling involves a weighted Delaunay triangulation of particle centers to form a tetrahedral mesh (Fig. 1). The solid volume within each tetrahedron is defined by the intersection of each tetrahedron with its vertex DEM spheres. Meanwhile, the fluid fraction consumes the remainder of each tetrahedron to form individual pores. Each pore is connected to four neighboring pores to constitute a pore network where a Stokes-flow is established based on an integral form of the continuity equation:

Θ i ∂ ρ f ∂t dV = - Θ i ∇ • (ρ f u)dV (3)
where Θ i is the domain of pore i, ρ f is the fluid density, and u is the fluid velocity. Application of the divergence theorem reduces the volume integral to a contour integral:

Θ i ∂ ρ f ∂t dV = - ∂ Θ i ρ f (u -v) • ndS - ∂ Θ i ρ f v • ndS (4)
where ∂ Θ i is the pore contour, v is the contour velocity, and n is the outward pointing unit vector. The consideration of fluid compressibility follows the relation of fluid bulk modulus, K, to the change of fluid pressure with respect to density:

K = ρ f ∂ P i ∂ ρ f (5)
where K can be a function of fluid pressure and air fraction as highlighted in Eq. 13. Finally, reducing ∂ Θ i to only the fluid fractions, S f i j , of the pore contour and assuming small Mach numbers, Eq. 4 becomes:

Θ i 1 K ∂ p i ∂t dV = 4 ∑ j=1 S f i j (u -v) • n dS -Vi (6)
where Vi is the rate of pore volume change and the integral on the right hand side represents the sum of fluid fluxes exchanged by each pore and its four neighbors ( j=1 to 4):

4 ∑ j=1 S f i j (u -v) • n dS = 4 ∑ j=1 q i j (7)
At low Reynolds numbers, the flux q i j through the pore throat connecting pore i and j is proportional to a local pressure gradient, with a coefficient k i j reflecting the local conductivity g i j :

q i j = g i j (p i -p j ) (8) Hence, 4 ∑ j=1 g i j (p i -p j ) = Vi + ṗi V i K (9)
which is the discrete form of a diffusion equation. Within Eq. 9, the pressures in the neighboring pores are represented by p i and p j and the length of the pore throat is l i j . The conductivity, g i j , was defined in [START_REF] Chareyre | Pore-Scale Modeling of Viscous Flow and Induced Forces in Dense Sphere Packings[END_REF]:

g i j = k d S f i j µ(T )l i j (10)
where k d is the permeability between pores i and j, S f i j is the area of the facet shared by pore i and j, and µ(T ) is the temperature dependent fluid viscosity.

The matrix representation of the full linear system representing Eq. 9 for all pores can be expressed as:

Gp = V ( 11
)
where G is the conductivity matrix, p is the vector containing the pressure within each pore, and V is the vector of rate of volume changes, which is composed of the volume changes due to: particle movements ( VG ), particle thermal expansion ( Vtm ), and fluid thermal expansion( Vth ):

V = VG + Vtm + Vth (12)
VG depends linearly on the particle velocities, ẋ, which can be expressed such that VG = Eẋ, where E is essentially a tensor with projected surface area.

In the case of gas-liquid mixtures, an equivalent fluid compressibility C eq is defined by following [START_REF] Najari | Thermo-hydro-mechanical response of granite to temperature changes[END_REF]:

C eq = φC a + (1 -φ )C w (13) C a = 1 P a (14) φ = P a,0 P a φ 0 ( 15
)
where φ is the fraction of gas within the mixture, C a is the compressibility of gas, P a is the absolute pressure of the gas and P a,0 and φ 0 are the initial absolute gas pressure and gas fraction respectively. C eq is updated at each time-step and used within the compressible flow scheme, Eq. 9, K = 1/C eq .

Heat transfer

Advection is simulated by treating each pore as an open thermodynamic system assuming that:

• radiation can be neglected (T<700 K)

• fluid kinetic energy is negligible (0<Re<1000)

• fluid compression does not generate heat Thus the integral form of the first law of thermodynamics becomes:

∆U o = Q + n ∑ m=1 ∆U m (16)
where U o is the internal energy of the open system (one pore in the present pore network), Q is a source term representing the quantity of heat supplied to the system, and ∆U m is the change of internal energy of the mth connected system (in the present pore network, these are neighboring pores and neighboring particles, Fig. 1). Considering the assumptions above, ∑ n m=1 ∆U m is limited to mass energy transfer, µ, and boundary heat flux, Φ. Thus, ∆U o can be expressed as a summation of surface integrals:

∆U o = Q + ∂ Θ i Φ • n dS + ∂ Θ i µ • n dS ( 17
)
where n is the outward pointing unit vector.

Advective heat transfer

Starting with the mass energy flux integral, each pore of the present pore network abuts four neighboring pores ( j=1 to 4), which means the integral can be reduced to a summation:

S µ • n dS = 4 ∑ j=1 S f i j µ i j • n dS = 4 ∑ j=1 µ i j ( 18
)
where the pore-pore mass-energy-flux between home pore i and neighbor pore j, µ i j , is the average over the shared facet and depends on the volumetric flow q i j (see Sec. 2.2):

µ i j = q i j c f ρ f T u (19) T u =      T i if q i j > 0 T j if q i j < 0      (20)
where T u is the temperature of the home or neighbor pore depending on the flow direction, c f , is the fluid heat capacity at constant pressure, and ρ f is the fluid density.

As presented in Sec. 2.2, the fluid flux q i j , is already solved (thus, the flow direction is known apriori)

as described in Sec. 2.2 as Gp. Therefore, the change of internal energy due to mass-energy-flux for all pores can be computed as:

U t+∆t i -U t i ∆t = c f ρ f 4 ∑ j=1 (T t u g i j ) p t+∆t i -p t+∆t j (21)
which, in matrix form, can be written as:

∆U = HP (22)
with

H = c f ρ f ∆t 4 ∑ j=1 (T t u g i j ) (23)
The internal energy, U t+∆t , is computed for all pores:

U t+∆t = U t + Hp (24)
and the temperature of the pores at t + ∆t is updated as:

T t+∆t = U t+∆t c f ρ f V t+∆t .
(25)

Conductive heat transfer

Conductive heat transfer is simulated between interacting particles, between neighboring pores, and between pores and particles assuming:

• radiation is neglected (T<700 K)

• the resistance to heat transfer inside the particle is significantly smaller than between particles (Biot

number = h i d i k i 1)
Given these assumptions, the heat conduction equation for a single particle k follows:

m k c p ∂ T k ∂t = ∇ • Φ k + Q (26)
where m k is the mass of particle k, c p is the particle heat capacity at constant pressure, T k is the temperature of the particle, Φ k is the boundary heat-flux into the particle and Q is the heat source supplied to the system. ∇ • Φ k is equivalent to the volume integral of the heat-flux, which can be reduced to a surface integral using the divergence theorem:

∇ • Φ k = ∂ Γ Φ k • ndS (27)
where ∂ Γ is the contour of particle k and n is the outward pointing unit vector. The surface integral is reduced to a summation along all pores and particles incident to particle k:

m k c p ∂ T k ∂t = N ∑ w=1 S p wk Φ p,wk • n wk dS + M ∑ u=1 A uk Φ s,uk • n uk dS (28)
where S p wk and Φ p,wk are the contact interface and average heat-flux, respectively, between particle w and k. and M and N are the number of pores and particles interacting with particle k, respectively.

Conductive heat-flux between contacting particles, Φ p,lm , follows existing methods by assuming the heat flux is linearly related to the inter particle temperature gradient by a thermal resistance (η) [START_REF] Feng | Discrete thermal element modelling of heat conduction in particle systems: Basic formulations[END_REF][START_REF] Liang | A new model for heat transfer through the contact network of randomly packed granular material[END_REF]:

S p lm Φ p,lm • n lm dS = η(T l -T m ) (29)
where T l is a unique temperature value for particle l and can be interpreted as the average particle temperature. η can be defined in various ways depending on the application of interest, as outlined in Sec. 3 and Sec. 4.

The first integral of Eq. 17, representing conductive heat-flux into the pore, is reduced to a summation by applying the divergence theorem:

∂ Θ i ΦdS = 4 ∑ j=1 S f i j Φ f ,i j • n i j dS + 4 ∑ k=1 A ik Φ s,ik • n ik dS (30)
where Φ f is the conductive heat flux within the fluid phase, n i j is the unit vector connecting pores i and j.

Φ s is the conductive heat flux between solid and fluid phases with n ik being the unit vector connecting pore i center to particle k. The second term of Eq. 30, representing the heat-flux between particles and pores (Φ s,ik ), matches the second term of the conservation of energy for each particle in Eq. 28. Therefore, the estimate is made once for each particle pore pair and follows a traditional spherical heat transfer approach [START_REF] Incropera | Fundamentals of Heat and Mass Transfer[END_REF], reducing the surface integral as follows:

A ik Φ s,ik • n ik dS = h ik A ik (T k -T i ) (31)
where A ik is the surface of particle k interacting with the pore i (spherical triangle shown in Fig. 1), T i is the temperature of the pore k, T k is the temperature of particle k, and h ik the heat transfer coefficient. h ik , is computed based on the Nusselt number Nu which can be empirically estimated using the macroscopic porosity of the particle assembly (0.35< ε <1) and Reynolds number (0<Re<10 2 ) as proposed by [START_REF] Tavassoli | Direct numerical simulation of fluid-particle heat transfer in fixed random arrays of non-spherical particles[END_REF]:

Nu k = (7 -10ε + 5ε 2 )(1 + 0.1Re 0.2 k Pr 1/3 k ) + (1.33 -2.19ε + 1.15ε 2 )Re 0.7 k Pr 1/3 k (32)
with Pr Prandtl's number and Re the Reynolds number based on the volume average of the pore k fluid velocity. The heat transfer coefficient then becomes

h ik = Nu k • k f /(2r i ),
with k f as the thermal conductivity of the fluid.

Given the particle-fluid (Φ s,ik , Eq. 31) and particle-particle (Φ p,i j , Eq. 59) heat flux approximations, Eq. 28 is approximated using a forward Euler scheme to estimate the particle temperature change:

T t+∆t k = ∆t m i c p N ∑ w=1 2r 2 c (k w + k k ) d wk (T t w -T t k ) + M ∑ u=1 Nu k k f A k 2r i (T t u -T t k ) + T t k ( 33
)
where M is the number of incident pores and N is the number of contacting particles.

Similarly to the particle conductive heat transfer, pores conduct heat between one another through their interface. Each interface is characterized by a thermal resistance, which determines the heat flux based on the temperature difference between pores, reducing the surface integral in Eq. 30 to:

S f i j Φ f ,i j • n i j dS = k f S f i j l i j (T j -T i ) (34) (35)
where l i j is the distance between the centers of pores i and j, T i and T j are the temperatures of the respective pores, and S f i j is the fluid area shared by both pores (Fig. 1).

In the present implementation, conservation of energy is ensured by computing pore-pore conduction using mid-step pore temperatures based on the advective heat flux. In other words, pore temperatures are updated twice per time step, first by advection, then by conduction.

Thermo-Hydro-Mechanical coupling

Hydro-Mechanical coupling

The two way hydro-mechanical coupling is only briefly described here -a full description and validation can be found in [START_REF] Chareyre | Pore-Scale Modeling of Viscous Flow and Induced Forces in Dense Sphere Packings[END_REF] and [START_REF] Scholtès | Modeling wave-induced pore pressure and effective stress in a granular seabed[END_REF]. The hydraulic force exerted by the pore fluid on the particles is decomposed into two contour integrals of the hydrostatic pressure P and viscous shear stress τ:

F p = Γ p P • ndS + Γ p τ • ndS = F P,p + F τ,p (36)
where Γ p is the solid surface of the particle p. Hydrostatic pressure forces, F P,p , and viscous forces F τ,p are both added to the interparticle force estimates (Eq. 1). Finally, the effect of particles movement on the pore network morphology is taken into account through Eq. 12.

Thermo-mechanical coupling

The thermo-mechanical coupling results from thermal expansion/contraction of both fluid and particles.

The associated pore volume changes are added to the existing rate of volume changes (Eq. 12) used within the flow solution (Eq. 11). In this way, the compressibility effects are handled inherently within the flow solver. The pore volume change contribution by the expansion/contraction of the solid phase ( Vtm ) depends entirely on solid temperature change (Sec. 2.3.2) and is expressed as:

∆r t+∆t k = r t k β p (T t+∆t k -T t k ) (37)
where β p is the coefficient of thermal expansion of the particles, r k is the radius of particle k, and T k its temperature. The particle volume changes, ∑ 4 k=1 ∆V i , are used to compute Vtm,i :

∆V i = - 4 ∑ k=1 A ik A k 4 3 π(r 3 k,t+∆t -r 3 k,t ) (38)
Vtm,i = ∆V i ∆t ( 39
)
where ∆t is the thermal time step, A ik is the spherical triangle area of particle k shared by pore i (Fig. ¿1), and and A k is the total surface area of particle k.

Meanwhile, the pore fluid thermal expansion/contraction is computed as:

Vth,i = V i β f (T )∆T ∆t ( 40
)
where ∆T is the temperature change and β f (T ) is the temperature dependent fluid volumetric coefficient of thermal expansion. Finally, Vth and Vtm contribute to the solution of fluid pressures (Eq. 12).

Porosity scaling

Various attributes can be scaled if the desired porosity is not equivalent to the exact packing porosity. For example, heat capacity, c p , is scaled such that that the total heat storage is representative of the desired material:

c p = c p,o 1 -φ d 1 -φ p (41)
where c p,o is non-scaled heat capacity of the desired material, φ d is the desired porosity, and φ p is the DEM sphere packing porosity.

Similarly, the pore space can be scaled such that such that the volume of fluid, V k , per tetrahedral matches the desired porosity:

V k = V k,i φ d φ p ( 42
)
where V k,i is the geometric pore volume, φ p is the porosity of the DEM sphere packing, φ d is the desired porosity.

Analytical verification of particle-particle conduction scheme

The numerical conduction scheme for particle-particle heat conduction (Sec 2.3.2) was verified by comparing various thermal DEM simulations (without pore fluid and neglecting mechanical interactions between particles) to the 1D heat equation:

∂ T ∂t = α ∂ 2 T ∂ x 2 (43)
The analytical solution to the 1D heat equation is constrained by the following initial and boundary conditions:

T (x, 0) = 120 x ∈ [0, L] (44)
T (0,t) = T (L,t) = 0 t ∈ [t > 0] (45)
Using Fourier series, the unsteady solution becomes:

T (x,t) = ∞ ∑ n=1 D n sin nπx L exp α n 2 π 2 t L 2 (46) D n = 2 L L 0 120 sin nπx L dx ( 47
)
where L ≈ 1 is the domain length, α is the effective thermal diffusivity, and t is the time.

Within the DEM conduction scheme, thermal resistivity estimate (η) emulates a continuum by modeling heat flux through a wall with depth d and area A depending on the interacting particles radii:

S p lm Φ p,lm • n lm dS = η(T l -T m ) (48) η = (k l + k m )/2 d A (49)
d = r l + r m (50) A = 4r l r m (51) (52)
Thermal micro-parameters of the DEM model are listed in Table 1. Readers can also find the practical input script for conduction in a single row of spheres titled "conductionVerification.py" in the supplementary material included with this paper.

Single row of spherical particles

A single row of DEM particles comprised of 166 particles was compared to the analytical solution for the 1D heat equation. The analytical solution was computed using an effective thermal diffusivity computed by scaling the density of DEM spherical elements to cubical continuum elements (ρ p π/6):

α = 6k p πC p ρ p . ( 53
)
The boundary particles were set to a constant temperature (T=0 • C) while the remaining particles were set to an initial temperature of 120 • C. The results show that the numerical temperature distributions match the analytical solution with space and time as shown in Fig. 2, with a maximum RMSE of 5.4e-06 for all Fourier numbers.

Random packings of spherical particles

The thermal diffusivity of random sphere packings cannot be estimated analytically. It can however be estimated numerically by simulating the 1D cooling scenario described in Sec. 3.1. Specifically, the thermal diffusivity α of the assembly is estimated considering the temperature evolution of three particles located along the axis of the packing (Fig. 3) as:

α = dT /dt d 2 T /dx 2 (54)
where the first and second derivatives are approximated by:

dT dt ≈ T t x -T t-∆t x ∆t (55) d 2 T dx 2 ≈ T t x+∆x -2T t x + T t x-∆x ∆x 2 (56)
The effective thermal diffusivity is monitored while the packing cools from an initial temperature of diffusivity for random packings was investigated by running 100 realizations of the thermal diffusivity test as shown in Figures 4 and5. Each packing realization contained the same uniform particle size distribution (ca. 970 particles, unif(0.003m±5%) or ca. 880 unif(0.003±30%)) but particle positions were randomly generated prior to each thermal diffusivity test. Following the 100 realizations, an estimated mean and standard deviation were assigned to the effective thermal diffusivity distribution of the material. As shown in Figures 4 and5, increasing the range of particle sizes corresponds to an increase of effective thermal diffusivity spread. These parameter estimates also enabled an illustration of temperature distribution variance with space and time (Figure 6). As shown, the variance is greatest where and when the heat flux is highest (i.e. at center of the specimen during the middle of the simulation). In summary, using random sphere packings reproduce analytical solutions within an expected statistical variation that depends on the sphere size distribution. Further, the thermal diffusivity estimate process presented is also the calibration process of the particle thermal conductivity according to the desired macroscopic thermal diffusivity. 2).

For both scenarios, the center body and center pore transient temperatures were compared. The THM-DEM and CFD models used identical material parameter values reported in Table 2. Readers can also find the practical input scripts for this section, titled "flowScenario.py" and "noFlowScenario.py", in Appendix C as well as the supplementary material included with this paper.

THM-DEM Thermal Resistivity

The conductive heat flux between two contacting DEM particles, Φ p,lm , was defined using the contact area, as presented in [START_REF] Norouzi | Coupled CFD-DEM Modeling[END_REF]:

Φ p,lm • n lm dS = η(T i -T m ) (57) η = 2r 2 c (k l + k m ) d lm (58)
d lm = r l + r m -p d (59) r c = 4d 2 lm r 2 m -(d 2 lm -r 2 l + r 2 m ) 2 2d lm (60) ( 61
)
where d is the distance between particles l and m less the overlap, p d , r is the particle contact radius, k is the particle thermal conductivity, and T l/m are the particle l and m temperatures.

Fully resolved solution

The fully resolved thermo-hydraulic CFD model was implemented in the commercial software package ANSYS CFX (CFX, 2019) which uses the Finite Volume Method to solve the governing equations presented below. Assuming laminar fluid flow (Re<2300), incompressible fluid, and no buoyancy forces the continuity equation follows:

∂ ρ ∂t + ∇ • (ρv) = S m ( 62
)
where ρ is density, v is velocity, t is time and S m is the mass source (e.g. due to vaporization of liquid droplets). The internal mass source was defined as null in the current model. Conservation of momentum in an inertial (non-accelerating) reference frame is described by

∂ ∂t (ρv) + ∇ • (ρvv) = -∇p + ∇ • ( τ) + ρg + F ( 63
)
where p is the static pressure and τ is the stress tensor. The gravitational body force, ρg and external body force, F were both neglected for the present comparison. The stress tensor τ is given by:

τ µ (∇ • v + ∇ • v T ) - 2 3 ∇ • vI ( 64
)
where µ is the molecular viscosity, I is the unit tensor, and the second term on the right-hand side is the effect of volume dilation.

The conservation of energy equation follows a low-speed flow variant:

∂ (ρe) ∂t + ∇ • (ρve) = ∇ • (λ ∇T ) -p∇ • v + ∇ • ( τ • v) + S h ( 65
)
where e is the internal energy and S h is an internal energy source. Due to very low fluid velocity and small pressure changes, terms p∇ • v and ∇ • ( τ • v) (viscous dissipation) were neglected as well as internal energy source S h . In consequence, energy change in time depends only on advection and conduction.

The transport equations were augmented with constitutive equations of state for density and enthalpy to form a closed system. It was assumed that fluid is incompressible and there are no buoyancy forces in the fluid. Hence, specific heat (at constant pressure) c p and density ρ are constant and the incompressible equation of state can be written

dH = c f dT + d p ρ ( 66
)
where h is enthalpy.

The coupling of pressure and velocity follows the high-resolution scheme based on discretization methods presented by [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF] and modified by [START_REF] Majumdar | Role of underrelaxation in momentum interpolation for calculation of flow with nonstaggered grids[END_REF].

Comparison of results

Fig. 9 shows the final comparison of body average temperatures at location (0.25,0.25,0.25) m. For the flow scenario, the THM-DEM model transfers heat more quickly to the center body and center pore which is due to the coarse fluid discretization of pores and particles. Meanwhile, the CFD model captures intricate movements of fluid around particles and complex heat transfer gradients within particles (Fig. 10).

Although the THM-DEM model is undoubtedly less accurate during loading (0-15 seconds), both solutions tend toward the same steady state solution. Computational comparison shows that THM-DEM is 100X faster than the ANSYS CFX solution. (Table 3). 2014)'s granite specimen (Fig. 12).

DEM Contact Model

The DEM contact model (Eq. 2, Sec. 2.1) follows a linear elastic model. The normal force, f n , between two interacting particles i and j is evaluated according to:

f n i j = k n i j ∆D i j • n n i j ( 67
)
where k n i j is the normal stiffness, n n i j is the unit vector parallel to the branch vector joining the centers of i and j, and ∆D i j = D i j -D eq i j is the displacement between i and j computed from the equilibrium distance D eq and the actual distance D.

k n is computed assuming two springs are in serial with lengths equal to the interacting particle radii:

k n = E a R a E b R b E a R a + E b R b ( 68
)
where E is a calibrated particle microparameter referred to as "micro Young's modulus" and R is the radius of particles a and b.

Since the shear force, f s , depends on the orientation of both particles, it is updated incrementally in a local coordinate system according to:

f s i j = f s i j,prev + ∆f s i j (69)
with

∆f s i j = k s i j ∆u s i j • n s i j ( 70
)
where k s i j is the shear stiffness, n s i j is the unit vector perpendicular to the branch vector, and ∆u s is the incremental tangential displacement.

DEM parameter calibration

The micro DEM properties reported in Table 4 were calibrated using typical DEM compression tests and the thermal diffusivity test discussed in Sec. 3. The thermal expansivity of the particles was calibrated by matching the macroscopic specimen thermal expansion to the experimentally observed thermal expansion for a given temperature change. Finally, a density scaling was used to increase the time-step of the explicit motion integration and thereby reduce the computational time to solution [START_REF] Itasca | PFC 2D-user manual[END_REF][START_REF] O'sullivan | Particulate discrete element modelling[END_REF][START_REF] Sheng | Numerical studies of uniaxial powder compaction process by 3D DEM[END_REF][START_REF] Thornton | Quasi-static shear deformation of a soft particle system[END_REF]. In the current quasi-static thermal-hydraulicmechanical heat transfer simulation, the propagation of elastic waves is the result of temperature changes Application of heat flux to particle packing yielding left) temperature and force on particle within packing and right) velocity of particle within packing for various particle densities only (expansion/contraction of DEM particles). Thus, it is ensured that the density scaling does not affect the development of particle forces. In the context of the present work, a density scaling analysis was performed by applying a ∆T =20 • C to one end of the specimen (Fig. 12) for 25 s while monitoring particle temperature, force (magnitude), and velocity. As shown in Fig 11, the density scaling affects the particle velocity, without affecting particle force. Thus, in the quasi-static problem presented here, a density scaling of 10 40 with a time step of 0.005 s remains stable without affecting the mechanical behavior of the system.

Hydro-Mechanical verification

The first experiment focused on deriving a permeability estimate based on the HM response of the granite specimen by pumping 3.333e-10 m 3 /s of water into the fluid filled cavity. Cavity pressure was monitored during cavity pressurization until steady state was achieved. After one hour at steady state, flow was terminated and depressurization within the cavity was monitored. As shown in Fig. 13, air fractions (φ 0 ) between 2.5e-5 and 3.75e-5 yielded the best match to experimental data. Similar to [START_REF] Najari | Thermo-hydro-mechanical response of granite to temperature changes[END_REF]'s FEM-HM model, higher φ 0 yields better accuracy during pressurization, and the lower φ 0 yields better accuracy during depressurization. Unlike [START_REF] Najari | Thermo-hydro-mechanical response of granite to temperature changes[END_REF], the air-fraction magnitude is roughly 2 orders of magnitude lower for the HM-DEM model. In both models, φ is a calibration parameter that accounts for more underlying complex physical interactions, so it is expected that the parameter does not match.

Thermo-Hydro-Mechanical verification

The second experiment presented by [START_REF] Najari | Thermo-hydro-mechanical response of granite to temperature changes[END_REF] focused on elucidating the THM response of the same granite specimen. During this experiment, the fully saturated granite rock specimen's surface was heated from 25 to 70 • C over the course of one hour (Fig. 14a). The temperature was held at 70

• C for three hours, after which the surface temperature was reduced to 25 • C and held for eight additional hours. During the course of this 12 hour thermal response test, the cavity pressure and temperature were monitored. As shown in Fig. 13, application of the same experimental boundary temperatures to the DEM-THM results in the same experimental and numerical pressurization and depressurization trends reported by [START_REF] Najari | Thermo-hydro-mechanical response of granite to temperature changes[END_REF]. Specifically, cavity pressure increases with increasing temperature followed by a decrease to 0 kPa once the maximum temperature is reached and steady-state is achieved. As soon as the temperature begins to drop, the cavity pressure decreases below 0 kPa and follows the reverse trend observed during heating. Numerical cavity temperatures also follow closely with experimental and numerical temperatures reported by [START_REF] Najari | Thermo-hydro-mechanical response of granite to temperature changes[END_REF], as shown in Fig. 14b.

k solid W m -1 K -1 30 3.0 3.0 C p J kg -1 K -1 1454 790 790 ρ kg m -3 2750e40 ca. 2750 2750 α m 2 s -1 - 1.43e-6 1.433e-6 β particle K -1
1.1e-4 3.0e-5 3.0e-5

β f luid K -1 β (T) - β (T) φ 0 2.5e-5 - 3e-3 E GPa

A Source code and installation

Figure 1 .

 1 Figure 1. Heat transfer model notations and geometric considerations for fluid flow, advection, and conduction models.

Figure 2 .

 2 Figure 2. Evolution of temperature distribution for a single row of DEM particles compared to the 1D analytical solution, with Fo h = αt L 2 .

Figure 3 .

 3 Figure 3. Temperature distribution within a random packing subjected to cooling at Fo h =0.06. Effective thermal diffusivity is estimated considering the temperature evolution at three points along the axis: T x+∆x , T x , T x-∆x

Figure 4 .

 4 Figure 4. 100 realizations of effective thermal diffusivity estimation associated for random packings generated with two different particle size distributions.

Figure 5 .Figure 6 .

 56 Figure 5. Distribution of effective thermal diffusivity estimate for two random packings generated with two different particle size distributions over 100 realizations.

Figure 7 .

 7 Figure 7. Comparison of left) DEM Yade sphere packing based pore network with right) ANSYS CFX mesh

Figure 8 .

 8 Figure 8. Specimen dimensions and boundary conditions for Yade DEM and ANSYS CFX a)no-flow scenario and b) constant-flow scenario.

Figure 9 .

 9 Figure 9. Temperature comparison for ANSYS CFX and Yade DEM in left) flow condition and right) no flow condition for body located at (0.024,0.028,0.026) and pore center located at (0.024,0.023,0.02545)

Figure 10 .

 10 Figure 10. Cross sectional temperature distribution for left) Yade DEM and right) ANSYS CFX at t=30s

  Figure11. Application of heat flux to particle packing yielding left) temperature and force on particle within packing and right) velocity of particle within packing for various particle densities

Figure 12 .

 12 Figure 12. Yade DEM specimen left) Isometric b) split to show fluid filled inner cavity

Figure 13 .

 13 Figure 13. Experimental and numerical cavity pressure curves a) permeability test b) thermal response test

  model for the simulation of heat transfer and associated thermo-hydro-mechanical couplings in particulate systems has been developed within the open source software Yade DEM. The proposed THM scheme enhances the capabilities of an existing and proven HM scheme by simulating the following heat transfer mechanisms: • conduction between contacting solid particles • conduction between pores • conduction between solid particles and pore fluid • advection through pore fluid flow Each mechanism can be simulated independently or collectively, as demonstrated through dedicated verification exercises. In particular, the successful comparison of the present THM-DEM model against a fully resolved CFD model proves that the geometrical considerations and heat flux models for the advection and conduction schemes are physically realistic despite their coarse treatment of the two phase domain. In addition to the verified conduction and advection schemes presented, thermo-mechanical couplings arise from temperature induced volumetric changes in both solid and fluid phases of the THM-DEM model.The thermo-mechanical effect is exhibited by verifying the full THM-DEM against a THM experiment involving the heating of a saturated rock sample. The proposed model reproduced the key internal fluid pressure trends observed during the THM experiment. First, fluid pressures increase within the rock material during surface temperature heating. Next, the pressure falls as steady state is achieved. Finally, the surface temperature is cooled and the fluid pressure decreases below initial pressures while steady-state is regained. Since the scheme is based on a triangulation of sphere packing poral space, the model applies more specifically to denser particle assemblies, such as those encountered in geomechanics. Thus, it opens up new possibilities for describing and understanding THM processes in porous media from a micromechanical viewpoint such as thermally induced microcracking in geothermal systems or transport process in rock and fractured rock masses.Zhou, Z. Y.,Yu, A. B., and Zulli, P. (2009). Particle scale study of heat transfer in packed and bubbling fluidized beds. AIChE Journal, 55(4):868-884. Šmilauer V. et al. (2015). Yade documentation 2nd ed. the yade project. Transport in Porous Media.

B#

  Setting up a Thermo-Hydro-Mechanical simulation in Yade Find the full list of available parameters and their full descriptions at https://yade-dem.org/doc/ yade.wrapper.html#yade.wrapper.ThermalEngine. Otherwise, a standard ThermalEngine simulation in Yade should be setup as follows: # add FlowEngine and ThermalEngine to typical engine list: O.engines=[ ... , FlowEngine(label="flow"), ThermalEngine(label="thermal"), VTKRecorder(recorders=[..., 'thermal', ... ]), # enable cavity model thermal.particleT0 = 20 thermal.particleK = 30. thermal.particleCp = heatCap thermal.particleAlpha = 3.0e-6 # solid expansion coeff thermal.particleDensity = density thermal.tsSafetyFactor = 0 thermal.uniformReynolds =10 # set a uniform reynolds number (only for very low vel.) → thermal.porosityFactor = 0.006/utils.porosity() thermal.tempDependentFluidBeta = True # fluid expansion coeff # impose cavity at desired locations: flow.imposeCavity((x,y,z)) # check temperature at desired locations: flow.getPoreTemperature((x,y,z)) # visualize fluid thermal quantities (temp, RE, etc.) flow.saveVtk('VTK/') *********************************************************************This program is free software; it is licensed under the terms of the * → # GNU General Public License v2 or later. See file LICENSE for details. * → # ************************************************************************* / # # Script demonstrating the use of ThermalEngine by comparing conduction → # scheme to analytical solution to Fourier (rod cooling with constant → # boundary conditions). See details in: # # Caulk, R., Scholtes, L., Kraczek, M., Chareyre, B. A # pore-scale Thermo-Hydro-Mechanical coupled model for particulate systems. #W/(mK) heatCap = 710. #J(kg K) t0 = 400. #K r = rad k = 2 * 2.0 * r # 2 * k * r Cp = 710. #triax, VTKRecorder(iterPeriod=500,fileName='VTK'+timeStr+identifier+'/spheres newton ] for b in O.bodies: if isinstance(b.shape, Sphere): b.dynamic=False# we only need flow engine to detect boundaries, there is no flow computed flow.boundaryUseMaxMin=[0,0,0,0,0,0] flow.thermalEngine=True flow.bndCondIsTemperature=[1,1,0,0,0,0] flow.thermalEngine=True flow.thermalBndCondValue=[0,0,0,0,0,0]

525##O#

  ***********************************************************************This program is free software; it is licensed under the terms of the * → # GNU General Public License v2 or later. See file LICENSE for details. * → # ************************************************************************* / # # Script demonstrating the use of ThermalEngine by monitoring still fluid → # temperature changes in a sphere packing. # Also serves as a validation script for comparison # with ANSYS CFD. See details in # Caulk, R., Scholtes, L., Kraczek, M., Chareyre, B. A # pore-scale Thermo-Hydro-Mechanical coupled model for particulate systems. → # # note: warnings for inifiniteK and Reynolds numbers = nan for boundary → # cells in regular packings are expected. It does not interfere with the → 05,0.05,0.05) # corners of the initial packing → thermalCond = 2. #W/(mK) heatCap = 710. #J(kg K) t0 = 333.15 #K ,fileName='VTK'+timeStr+identifier+'/sphere newton ] for b in O.bodies: if isinstance(b.shape, Sphere): b.dynamic=False # mechanically static for b in O.bodies: if isinstance(b.shape, Sphere): dist = b.state.pos -Vector3(x,y,z) if np.linalg.norm(dist) < np.linalg.norm(cDist): bodies along x axis axis = np.linspace(mn[0], mx[0], num=5) axisBodies = [None] * len(axis) axisTrue = np.zeros(len(axis)) for i,x in enumerate(axis): axisBodies[i] = bodyByPos(x, mx[1]/2, mx[2]/2) axisTrue[i] = axisBodies[i].state.pos[0] from yade import plot Flow example script: 526 # ***********************************************************************This program is free software; it is licensed under the terms of the * → # GNU General Public License v2 or later. See file LICENSE for details. * → # ************************************************************************* / # # Script demonstrating the use of ThermalEngine by permeating warm fluid → # through a cold packing. Also serves as a validation script for comparison → # with ANSYS CFD. See details in # Caulk, R., Scholtes, L., Kraczek, M., Chareyre, B. A # pore-scale Thermo-Hydro-Mechanical coupled model for particulate systems. → # # note: warnings for inifiniteK and Reynolds numbers = nan for boundary → # cells in regular packings are expected. It does not interfere with the 05,0.05,0.05) # corners of the initial packing → thermalCond = 2. #W/(mK) heatCap = 710. #J(kg K) t0 = 333.15 #K

Table 1 .

 1 DEM Microparameters

	Parameter Value	units
	k p	2.0	W/(m K)
	C p	710	J/(kg K)
	ρ p	2600	kg/m 3
	r	0.003	m

Table 2 .

 2 DEM and CFD parameters

	Solid Parameter	Value
	k W/(m K)	2
	C p J/(kg K)	710
	ρ kg/m 3	4976
	T 0	• C	60
	Fluid parameter	Value
	k W/(m K)	0.65
	C f J/(kg K)	4184
	ρ kg/m 3	1000
	T 0	• C	60
	viscosity Pa • s	0.001
	Reynolds	0-530
	Incompressible
	Boundary conditions
	Fluid -Y dirichlet • C	70
	Fluid -Y dirichlet Pa	10
	Fluid +Y dirichlet Pa	0
	Fluid dirichlet ±X±y±z Pa	70
	Simulation duration (seconds)	30
	Body temp comparison location (0.25,0.25,0.25)
	Flux (kg/s) DEM	0.00394
	Flux (kg/s) CFD	0.00472
	begin

Table 3 .

 3 Performance comparison for THM-DEM and ANSYS CFX

		No-flow (hrs) Flow (hrs)	Cells
	Pore scale THM-DEM
	(Yade git-dc2ecaec (10-core))	0.04	0.06	1000
	Fully resolved FVM
	(ANSYS CFX v.19.2 (12-core))	3.78	5.38	3.5 million

Table 4 .

 4 DEM Calibrated microparameters and emergent macroparameters Parameter DEM micro DEM macro FEM/Experimental

  Readers are encouraged to contact the Yade community at https://answers.launchpad.net/ yade with questions regarding installation, usage, modifications, and theory.

	Source code is freely available as part of Yade DEM on gitlab.com. Installation on a Ubuntu linux
	(https://yade-dem.org/doc/installation.html) requires a single command:
	sudo apt-get install yade
	However, the source code for Yade DEM and the included thermal components presented here are available
	online for review and modification:
	• Full source code:https://gitlab.com/yade-dev/trunk
	• Thermal	component:https://gitlab.com/yade-dev/trunk/-/blob/master/
	pkg/pfv/Thermal.cpp

Acknowledgments

The material presented is based upon research supported by the Chateaubriand Fellowship of the Office for Science & Technology of the Embassy of France in the United States. Support was also provided by Veronica Eliasson. The fourth author would like to thank support from Tech21. The third author would like to thank support from the project "Fracture propagation in rocks during hydro-fracking -experiments and discrete element method coupled with fluid flow and heat transport" (years 2019-2022) financed by the National Science Centre (NCN) (UMO-2018/29/B/ST8/00255).