M.D Rafi Sonia
email: sonia.rafi@chu-rennes.fr

M.D Jean-Marc Tadie
email: jeanmarc.tadie@chu-rennes.fr

M.D Gagouin Arnaud
email: arnaud.gacouin@chu-rennes.fr

M.D Leurent Guillaume
email: guillaume.leurent@chu-rennes.fr

M.D Bedossa Marc
email: marc.bedossa@chu-rennes.fr

L E Tulzo

M.D Maamar Adel
email: adel.maamar@chu-rennes.fr

M.D Sonia Rafi

Doppler sonography of cerebral blood flow for early prognostication after out-of-hospital cardiac arrest: DOTAC study

Keywords: Out-of-hospital cardiac arrest, transcranial Doppler sonography, cerebrovascular circulation, prognosis

Aim: To assess the neurological prognosis of comatose survivors of cardiac arrest by early transcranial Doppler sonography (TCD). Methods: This was a prospective study performed between May 2016 and October 2017 in a medical intensive care unit (ICU) and a cardiac ICU of a university teaching hospital. All patients older than 18 years who were successfully resuscitated from an out-of-hospital cardiac arrest (OHCA) with persistent coma after the return of spontaneous circulation (ROSC) were eligible. We excluded patients for whom OHCA was associated with traumatic brain injury, no possibility of TCD measurements, or who were dead before establishing the neurological prognosis.

We measured the pulsatility index (PI) and diastolic flow velocity (DFV) of the right and left middle cerebral arteries within 12 h after ICU admission. The lowest DFV and highest PI values were used for the statistical analysis. The neurological outcome at hospital discharge was evaluated by the cerebral performance category. Results: Forty-two patients were included in the final analysis: 15 had good and 27 poor neurological outcomes. The PI was higher in the poor outcome (1.49 vs. 1.12, p = 0.01) than good outcome group and the DFV was lower in the poor outcome group (17.3 cm.s -1 vs. 26.0 cm.s -1 ; p = 0.01). Conclusion: Data provided by early TCD after ROSC are associated with neurological outcome. The use of TCD could help guide interventions to improve cerebral perfusion after ROSC in patients resuscitated from OHCA.

A C C E P T E D M A N U S C R I P T Introduction

Out-of-hospital cardiac arrest (OHCA) is a major cause of death and disability, with an incidence rate of 38 to 55 OHCA per 100,000 person-years worldwide [START_REF] Gräsner | EuReCa ONE27 Nations, ONE Europe, ONE Registry[END_REF]. Despite recent improvements in the management of OHCA, the survival rate at hospital discharge is still below 10% [START_REF] Bossaert | European Resuscitation Council Guidelines for Resuscitation 2015[END_REF]. Two types of events have been shown to be responsible for the death of patients who survive the initial phase of prehospital care: early reperfusion syndrome, which can induce a refractory state of shock, with death within the first 24 h following admission, and delayed death induced by neurological impairment [START_REF] Mongardon | Postcardiac arrest syndrome: from immediate resuscitation to long-term outcome[END_REF]. A retrospective analysis found that more than 60% of such patients died due to neurological injury [START_REF] Laver | Mode of death after admission to an intensive care unit following cardiac arrest[END_REF]. Neurological injury is not only the result of the anoxic period of cardiac arrest, but can also be exacerbated during the postresuscitation phase due to modifications in cerebral blood flow (CBF).

Transcranial Doppler sonography (TCD) is a non-invasive bedside technique that uses ultrasound to measure cerebrovascular haemodynamics [START_REF] Saqqur | Role of transcranial Doppler in neurocritical care[END_REF]. The transtemporal measurement of the middle cerebral artery (MCA) pulsatility index (PI) provides a surrogate of cerebral perfusion by recording the velocity of red blood cells circulating in this vessel [START_REF] Geeraerts | Le doppler transcrânien en neurochirurgie[END_REF]. A higher PI is induced by an increase in downstream vascular resistance [START_REF] Vigue | Le doppler transcrânien en réanimation[END_REF]. TCD is not currently part of the algorithm recommended in the neurological prognosis strategy for comatose survivors of cardiac arrest [START_REF] Sandroni | Prognostication in comatose survivors of cardiac arrest: An advisory statement from the European Resuscitation Council and the European Society of Intensive Care Medicine[END_REF].

Methods

Study design

This prospective observational study was conducted between May 2016 and October 2017 in the medical intensive care unit (ICU) and cardiac ICU of a university teaching hospital. The hospital's ethical committee approved this study (no. 16.34). The ethical committee waived informed consent, as TCD is considered to be observational.

Patient eligibility

All adult patients that were successfully resuscitated from an OHCA with persistent coma after the return of spontaneous circulation (ROSC) were eligible. Exclusion criteria were associated traumatic brain injury, the inability to perform TCD because of the absence of an acoustic window or because > 12 h had passed since admission, and death from refractory shock within the first 72 h, prior to neurological assessment.

OHCA management

All patients were managed upon admission based on recommendations and received continuous administration of intravenous midazolam and morphine with the objective of attaining a Richmond Assessment Sedation Scale (RASS) of minus five [START_REF] Bossaert | European Resuscitation Council Guidelines for Resuscitation 2015[END_REF][START_REF] Sessler | The Richmond Agitation-Sedation Scale: Validity and Reliability in Adult Intensive Care Unit Patients[END_REF]. Such sedation was maintained for 24 h or more if there was any reason for maintaining deep sedation, especially in cases of severe acute respiratory distress syndrome. All patients received targeted temperature management during the first 24 h with an objective of 36°C [START_REF] Nielsen | Targeted Temperature Management at 33°C versus 36°C after Cardiac Arrest[END_REF]. Arterial line placement for all patients allowed continuous monitoring of arterial pressure. A coronary Neuroprotective measures to prevent secondary cerebral injury of systemic origin were implemented [START_REF] Moeschler | Concept d'agression cérébrale secondaire d'origine systémique (ACSOS)[END_REF]. Neither department performed invasive monitoring of intracranial pressure.

Data collection and TCD measurements

The following parameters were recorded for each eligible patient upon admission: age, gender, mechanism of cardiac arrest (asystole, ventricular fibrillation, pulseless electrical activity), and the time interval between cardiac arrest and resuscitation manoeuvres (no flow), the start of resuscitation manoeuvres and ROSC (low flow), and ROSC and TCD examination.

TCD measurements were performed within 12 h after ICU admission using a 5-MHz probe for all patients (Philips Healthcare TM , CX50). We measured systolic flow velocity (SFV), diastolic flow velocity (DFV), mean flow velocity (MFV), and PI of the MCA of each patient for both sides. The PI is a parameter used to assess pulsatility and is defined as the difference between maximum and minimum flow velocity as follows: PI = (SFV-DFV)/MFV [START_REF] Geeraerts | Le doppler transcrânien en neurochirurgie[END_REF]. Two trained operators performed all measurements.

We recorded parameters that can influence TCD values during measurement [START_REF] Vigue | Le doppler transcrânien en réanimation[END_REF]: the arterial partial pressure of carbon dioxide (mmHg); mean, systolic, and diastolic blood pressure (mmHg); heart rate; core temperature (°C); haemoglobin level (g.dl -1 ); and treatment with vasopressors (i.e., epinephrine and/or norepinephrine) at any dose. Serum lactate levels and oxygen partial pressure during the TCD measurements were also measured. If one or more of the biological parameters were not recorded concurrently with the Doppler measurements, the results of the parameter(s) measured closest to the time of the TCD were collected.

A C C E P T E D

M A N U S C R I P T
The neurological outcome at hospital discharge was evaluated according to the Pittsburgh Cerebral Performance Category (CPC) scale [START_REF] Gehlbach | Nicotine absorption by workers harvesting green tobacco[END_REF][START_REF] Hazinski | Recommended guidelines for reviewing, reporting, and conducting research on in-hospital resuscitation: the in-hospital 'Utstein style[END_REF] by two independent intensivists blinded to the TCD values. If there was a discrepancy, a third intensivist was asked to decide and validate the final result.

Neurological outcomes were divided into good and poor. A good neurological outcome was defined as a CPC of 1 (good recovery) or 2 (moderate disability). A poor neurological outcome was defined as a CPC of 3 (severe disability), 4 (vegetative state), or 5 (death due to a comatose state or brain death) [START_REF] Mongardon | Postcardiac arrest syndrome: from immediate resuscitation to long-term outcome[END_REF][START_REF] Laver | Mode of death after admission to an intensive care unit following cardiac arrest[END_REF]. Patients who died before that neurological evaluation was possible were excluded from the analysis, as no neurological evaluation could be performed on sedated patients.

The N20 response to early-latency evoked potentials was sought in accordance with the recommended algorithm for determining the neurological prognosis of comatose survivors of a cardiac arrest [START_REF] Sandroni | Prognostication in comatose survivors of cardiac arrest: An advisory statement from the European Resuscitation Council and the European Society of Intensive Care Medicine[END_REF] and a decision to withdraw and withhold life-sustaining healthcare was collectively made after a multimodal evaluation according to this algorithm and blinded to the TCD results.

Statistical analysis

A minimum sample of 20 patients for each group was calculated to detect a difference of 0.5 between PI values of each group with a power of 0.8 and an α risk < 0.05 between patients with good and poor neurologic outcomes [START_REF] Pollock | Anoxic Injury-Associated Cerebral Hyperperfusion Identified with Arterial Spin-Labeled MR Imaging[END_REF].

Continuous variables are expressed as medians (interquartile range) and categorical variables as numbers (%). The Mann-Whitney U-test or Fisher's exact tests were applied, as appropriate, for comparisons between groups. Statistical analyses were performed using Statview 5.0 (SAS Institute Inc., Cary, NC). A p-value < 0.05 was considered significant.

A C C E P T E D

M A N U S C R I P T

Results

Study population

During the study period, 88 cardiac-arrest patients resuscitated with ROSC were screened, 36 with in-hospital cardiac arrest and 52 with OHCA (Figure 1). Among the 52 patients with OHCA, 10 (19.3%) were excluded (two because of the absence of an acoustic window at TCD examination, five because of death before neurological evaluation, and three because > 12 h had passed since admission). Finally, 42 patients (80.7%) met the eligibility criteria. The demographic characteristics of the study population are summarized in Table 1. Among these 42 patients, 15 (35.7%) had a good neurological outcome (11 with CPC = 1 and four with CPC = 2), and 27 (64.3%) had a poor neurological outcome (five with CPC = 3; 11 with CPC = 4; and 11 with CPC = 5) at discharge. The baseline characteristics, including demographics and haemodynamic parameters did not differ between the good and the poor neurological outcome groups (Table 2).

Poor neurological outcome

Eleven (26%) patients died before 72 h from neurological causes and were categorized in the poor outcome group. The N20 response was sought for five patients, including a bilateral absence for four. Those patients with a bilaterally absent N20 response died within the first 72 h after withdrawal of life-sustaining therapies. One had an intact N20 response but multiple comorbidities, including amyotrophic lateral sclerosis, so that the decision to withdraw lifesustaining therapies was made at 72 h. Among the six others, five were in a state of clinical brain death with persistent coma, the absence of brainstem reflexes, and the absence of a breathing drive during apnea, apart from any residual sedation test. The last had a CPC score of 4. Life-sustaining therapies were withdrawn at 72 h for these six patients.

A C C E P T E D

M A N U S C R I P T
Of note, five patients (9.6%) died prematurely from multiple organ failure while they were still sedated without neurological evaluation and consequently excluded from the final analysis.

These patients had TCD values within the normal range with a median IP of 1.16 (0.92-1.4) and median DFV of 30.3 cm.s -1 (18.9-41.7).

TCD measurements

TCD was performed after admission in the ICU. SFV, DFV, and MFV were recorded. The PI was automatically calculated by the ultrasound software as follows: (SFV-DFV)/MFV. The median time for TCD measurement after ROSC was 330 min (249-491) for the good-outcome group and 303 min (26-464) for the poor-outcome group (p = 0.93). The TCD values were significantly different depending on the neurological outcome (Table 2). The PI was higher in the poor-outcome (1.49 vs. 1.12, p = 0.01) than good-outcome group and the DFV was lower in the poor outcome group (17.3 cm.s -1 vs. 26.0 cm.s -1 ; p = 0.01).

Discussion

We found that patients with poor neurological outcomes after OHCA had a lower DFV and higher PI at admission than patients with good neurological outcomes, despite no differences in systemic arterial pressure.

Various studies have provided conflicting results concerning the association between initial

TCD values and neurological outcome. Wessels et al. [START_REF] Wessels | The prognostic value of early transcranial doppler ultrasound following cardiopulmonary resuscitation[END_REF] performed serial TCD examinations on various vessels (anterior, middle, and posterior cerebral arteries) during the first 72 h of hospitalization in 39 patients with cardiorespiratory arrest. Of the 39 patients, 22 died or had neurological outcomes considered to be unfavourable. By the fourth hour, the peak systolic and diastolic velocities in the middle cerebral artery were lower in non-survivors than survivors. However, patients with unfavourable neurological outcomes appeared to have higher resistance indexes, even if the result did not reach the significance. To date, only one study has shown a

A C C E P T E D M A N U S C R I P T
significant association between values measured by TCD examination and neurological evolution measured by the CPC score. Interestingly, the TCD examination was performed early in the course of ROSC, at three time points: before targeted temperature management, between 12 and 24 h, and between 24 and 48 h during the targeted temperature management phase. However, the study population consisted exclusively of children [START_REF] Lin | Transcranial Doppler ultrasound in therapeutic hypothermia for children after resuscitation[END_REF].

The main difference between these studies and ours is that we excluded patients who died prematurely from multi-organ failure without neurological evaluation. Of note, patients who died prematurely in our study had normal TCD values.

Concerning haemodynamic parameters, we found no differences in systolic, diastolic, or mean blood pressure (MBP) between the good-outcome and poor-outcome groups. In particular, only the MBP value did not differentiate patients whose cerebral perfusion was altered, contrary to the TCD parameters. We did not find a statistical difference in arterial blood pressure between our two groups, although we cannot rule out a clinically relevant difference due to the small sample.

The neurological outcome depends on the prompt restoration of systemic circulation and oxygenation after ROSC to meet cerebral oxygen demands [START_REF] Lin | Transcranial Doppler ultrasound in therapeutic hypothermia for children after resuscitation[END_REF][START_REF] Álvarez-Fernández | Use of transcranial Doppler ultrasound in the management of post-cardiac arrest syndrome[END_REF]. These studies suggest that CBF and the MCA pulsatility index at the early ROSC phase may correlate with neurological outcome. The high PI associated with decreased DFV observed in our study favours an increase in vascular resistance at the acute phase after a resuscitated cardiac arrest, as observed by other authors [START_REF] Buunk | Cerebral vasoconstriction in comatose patients resuscitated from a cardiac arrest?[END_REF][START_REF] Lemiale | Changes in cerebral blood flow and oxygen extraction during post-resuscitation syndrome[END_REF]. During the post-resuscitation phase after ROSC, ischaemia-reperfusion phenomena induce disturbances of CBF. Recent studies have described the course of DSC after ROSC, which can be divided into three periods [START_REF] Van Den Brule | Cerebral Perfusion and Cerebral Autoregulation after Cardiac Arrest[END_REF][START_REF] Iordanova | Alterations in Cerebral Blood Flow after Resuscitation from Cardiac Arrest[END_REF]. Immediately after ROSC, CBF is characterised by early hyperaemia for 30 min, followed by a hypoperfusion phase of 6 to 12 h. The third period shows restoration of normal CBF, increased CBF [START_REF] Pollock | Anoxic Injury-Associated Cerebral Hyperperfusion Identified with Arterial Spin-Labeled MR Imaging[END_REF][START_REF] Iida | Delayed_hyperemia_causing_intracranial_hypertension_after_cardiopulmonary.html[END_REF], or persistent

A C C E P T E D M A N U S C R I P T
decreased CBF from 12 to 72 h after ROSC [START_REF] Pollock | Anoxic Injury-Associated Cerebral Hyperperfusion Identified with Arterial Spin-Labeled MR Imaging[END_REF][START_REF] Iida | Delayed_hyperemia_causing_intracranial_hypertension_after_cardiopulmonary.html[END_REF][START_REF] Buunk | Cerebral vasoconstriction in comatose patients resuscitated from a cardiac arrest[END_REF]. There are several hypotheses that can explain such CBF abnormalities, including the secretion of endothelin, a potent vasoconstrictor that may be involved in the modification of cerebral vascular tone [START_REF] Bacic | Secretion of immunoreactive endothelin-1 by capillary and microvascular endothelium of human brain[END_REF]. Until recently, it was assumed that under normal circumstances, autoregulation maintains a constant CBF and changes in mean blood pressure within a range of 50-150 mm Hg have a minor influence on CBF [START_REF] Aaslid | Cerebral autoregulation dynamics in humans[END_REF][START_REF] Lassen | Cerebral Blood Flow and Oxygen Consumption in Man[END_REF]. However, recent studies have shown that autoregulation maintains CBF within a smaller range above baseline MBP [START_REF] Van Den Brule | Cerebral Perfusion and Cerebral Autoregulation after Cardiac Arrest[END_REF][START_REF] Tan | Defining the characteristic relationship between arterial pressure and cerebral flow[END_REF]. Several studies showed that cerebral autoregulation is preserved in the aftermath of cardiac arrest [START_REF] Buunk | Cerebral vasoconstriction in comatose patients resuscitated from a cardiac arrest?[END_REF][START_REF] Bisschops | Preserved metabolic coupling and cerebrovascular reactivity during mild hypothermia after cardiac arrest*[END_REF]. In contrast, several lines of evidence have demonstrated that there is no regulation of CBF. For example, Sundgreen et al. [START_REF] Sundgreen | Autoregulation of cerebral blood flow in patients resuscitated from cardiac arrest[END_REF] showed that cerebral autoregulation was either absent or altered, with the curve shifted to the right, in most of a series of 18 patients. Therefore, systemic arterial pressure may not allow prompt monitoring in brain perfusion. In the event of impairment of cerebral autoregulation, minimal variations in MBP are likely to alter CBF without being detectable by conventional monitoring, including continuous monitoring of intracranial pressure.

Moreover, it may be important to determine whether TCD values at an early stage after ROSC are predictive of neurological outcome [START_REF] Lemiale | Changes in cerebral blood flow and oxygen extraction during post-resuscitation syndrome[END_REF][START_REF] Doepp | Duplex sonography of cerebral blood flow after cardiac arrest-A prospective observational study[END_REF]. Impaired DFV and increased PI are considered to be pathological in severe traumatic brain injury [START_REF] Ract | Transcranial Doppler ultrasound goal-directed therapy for the early management of severe traumatic brain injury[END_REF] and are used to guide therapy to restore adequate cerebral perfusion. Therefore, normalization of TCD values at early stages of ROSC may represent an interesting approach to improve neurological outcome after OHCA.

Heradstveit et al. showed the presence of ischaemic lesions at 24 h in these patients by magnetic resonance imaging with diffusion sequences [START_REF] Heradstveit | Repeated magnetic resonance imaging and cerebral performance after cardiac arrest-A pilot study[END_REF]. The final objective would be to improve cerebral perfusion by optimizing local factors (parameters of mechanical ventilation, metabolic coupling, patient posture), as well as systemic factors (stroke volume, myocardial contractility).

Only a small number of patients were excluded because of the impossibility of interpreting the TCD (no window or the passage of too much time after admission). This confirms the feasibility

A C C E P T E D M A N U S C R I P T
of the TCD examination in common practice in ICUs during the acute phase after a cardiac arrest. The other exclusions consisted of five patients (9.6%) who died prematurely from multiple organ failure while they were still sedated and could not be neurologically evaluated. This finding is in agreement with the literature, which shows that approximately 10% of patients resuscitated from an OHCA early die of multiple organ failure [START_REF] Laver | Mode of death after admission to an intensive care unit following cardiac arrest[END_REF]. Such early death should not be considered for the neurological prognosis because the outcome could have been favourable if the patients had survived the initial acute circulatory failure.

Our study had some limitations. First, TCD does not directly measure CBF. Velocities may reflect CBF, provided there is no variation in the calibre of the insonated vessel. The relationship between TCD velocities and CBF are complex and involve many parameters, such as vasoreactivity to hypoxia and hypercapnia, the effect of neuronal metabolism, and the effect of the autonomic nervous system [START_REF] Willie | Integrative regulation of human brain blood flow: Integrative regulation of human brain blood flow[END_REF][START_REF] Hoedemaekers | Low cerebral blood flow after cardiac arrest is not associated with anaerobic cerebral metabolism[END_REF].

Second, the patients came from two different settings of the same hospital, the results of which merit confirmation on a larger scale. Nevertheless, the management of comatose survivor patients from an OHCA in the ICU is standardized due to the standard edition of international recommendations on the subject [START_REF] Nolan | European Resuscitation Council and European Society of Intensive Care Medicine Guidelines for Post-resuscitation Care 2015[END_REF]. Our population was small and heterogenous and a larger study is needed to perform multivariate analyses to identify PI or DVF as independent factors of a poor or good outcome. Furthermore, a larger sample would contribute to a better understanding of the pathophysiology of post-anoxic encephalopathy by providing data from patients with cardiac versus non-cardiac causes of OHCA.

Third, we acknowledge a potential bias given the low interrater reliability of some TCD measurements. All TCD measurements in our study were thus performed by only two trained operators to mitigate such bias.

A C C E P T E D

M A N U S C R I P T
Only one study has had a larger number of patients than ours; however, it did not show any relationship between TCD values and neurological outcome. This may be because the TCD measurements were performed later, i.e., at 48 h, 72 h, and 7 days after ROSC, including on patients who died before neurological evaluation [START_REF] Doepp | Duplex sonography of cerebral blood flow after cardiac arrest-A prospective observational study[END_REF].

Conclusion

Data provided by early TCD examination after ROSC are associated with neurological outcome. The use of TCD could help guide interventions to improve cerebral perfusion after ROSC in patients resuscitated from OHCA.

Conflicts of interest

  Based on the hypothesis that impaired CBF in the hours following an OHCA is a marker of ongoing cerebral injury, we assessed the relationship between CBF observed by early TCD examination and neurological outcome at hospital discharge.

  by percutaneous coronary intervention, was performed if the patient's attending physician deemed it necessary.

Figure 1 .

 1 Figure 1. Flow chart of patients study

Acknowledgements

The authors would like to thank the ICU staff and cardiology department for their participation in the transcranial Doppler measurements and their involvement in allowing the inclusion of patients admitted after OHCA.

Conflict of Interest Statement

All authors have made contribution to the conception of the study, acquisition of data, and analysis of data, and state the absence of conflict of interest in their contribution to this study


Martin Schmelzer

Richard Dwight

Paola Cinnella

Data-Driven Deterministic Symbolic Regression of Nonlinear Stress-Strain Relation for RANS Turbulence Modelling

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. Introduction

Turbulence modelling is a key-challenge for computational fluid dynamics (CFD) especially in industry. Despite considerable progress made in the field of high-fidelity turbulence modelling, such as large-eddy simulation (LES) and direct numerical simulation (DNS), RANS continues to be the standard approach used to predict a wide range of flows [START_REF] Slotnick | CFD Vision 2030 Study: A Path to Revolutionary Computational Aerosciences[END_REF]. However, using the less-computationally demanding RANS approach comes at the price of uncertainty due to approximate physical modelling of the turbulence closure.

Despite Reynolds-stress models (RSM), based on transport equations for the Reynolds-stress tensor components, offer the more elaborate second-moment closures for RANS, modelling of industrial flows is mainly done using linear eddy viscosity models (LEVM) [START_REF] Basara | A new hybrid turbulence modelling strategy for industrial CFD[END_REF][START_REF] Pope | Turblent Flows[END_REF][START_REF] Leschziner | Statistical Turbulence Modelling for Fluid Dynamics -Demystified: An Introductory Text for Graduate Engineering Students[END_REF]. This is due to the fact that RSM involve a more complex mathematical structure, which requires additional modelling and introduces computational challenges, whilst the method does not offer a superior predictive capability for all flows [START_REF] Leschziner | Statistical Turbulence Modelling for Fluid Dynamics -Demystified: An Introductory Text for Graduate Engineering Students[END_REF]. However, for flows with streamline curvature, adverse pressure gradients, flow separation or rotation LEVM do not deliver reliable predictions due to their inherent inability to predict the anisotropy of the Reynolds-stress correctly. Explicit Algebraic Reynolds-stress Models (EARSM), first introduced by Pope [START_REF] Pope | A more general effective-viscosity hypothesis[END_REF] and further developed by Gatski and Speziale [START_REF] Gatski | On explicit algebric stress models for complex turbulent flows[END_REF], have the potential to fill the gap by offering higher predictive fidelity compared to LEVM and being numerically more robust than RSM at similar computational costs as LEVM [START_REF] Wallin | Engineering turbulence modelling for CFD with a focus on explicit algebraic Reynolds stress models by[END_REF]. EARSM are derived from a projection of (simplified) RSM onto a set of tensorial polynomials (see below) [START_REF] Pope | Turblent Flows[END_REF][START_REF] Leschziner | Statistical Turbulence Modelling for Fluid Dynamics -Demystified: An Introductory Text for Graduate Engineering Students[END_REF], leading to a non-linear stress-strain relation for the Reynolds-stress tensor. As such, these models can be seen as higher-order extensions of the linear eddy viscosity concept. The main simplification is the omission of anisotropy transport in order to enhance the numerical robustness, which however also reduces the predictive potential of this modelling strategy [START_REF] Leschziner | Statistical Turbulence Modelling for Fluid Dynamics -Demystified: An Introductory Text for Graduate Engineering Students[END_REF].

Recently a new approach based on symbolic regression utilising genetic programming (GP) was introduced to learn the non-linear stress-strain relationship for the anisotropy tensor based on high-fidelity data from DNS and LES [START_REF] Weatheritt | A novel evolutionary algorithm applied to algebraic modifications of the RANS stress-strain relationship[END_REF][START_REF] Weatheritt | The development of algebraic stress models using a novel evolutionary algorithm[END_REF]. This data-driven method retains the input quantities used to derive EARSM but replaces the commonly used projection method to find the formal structure of the model by an evolutionary process based on model fitness. In that way it produces models similar to EARSM but with a mathematical form proven to reproduce the data it was trained on. This method has the potential to generate numerically robust models with a high predictive fidelity. Even though the method is non-deterministic it discovers similar expressions for different runs. However, it is not clear if this variability comes from the data or is due to the inherent randomness of GP.

To overcome this characteristic of GP non-evolutionary methods for symbolic regression have been introduced, e.g. Fast Function Extraction (FFX) [START_REF] Mcconaghy | FFX: Fast, Scalable, Deterministic Symbolic Regression Technology[END_REF], Elite Bases Regression (EBR) [START_REF] Chen | Elite Bases Regression: A Real-time Algorithm for Symbolic Regression[END_REF], Sparse identification of nonlinear dynamics (SINDy) [START_REF] Brunton | Discovering governing equations from data: Sparse identification of nonlinear dynamical systems[END_REF] or PDE functional identification of nonlinear dynamics (PDE-FIND) [START_REF] Rudy | Data-driven discovery of partial differential equations[END_REF]. These methods are all based on sparsity-promoting linear regression and show for a couple of cases similar or better performance and higher convergence rates for high-dimensional problems than GP. Due to their deterministic nature, they discover always the same model given the input quantities and parameters. Additionally, by varying the input parameters of the method, a hierarchy of models of different complexity and predictive fidelity are discovered, which can be used to assess overfitting for prediction.

In this work a deterministic version of symbolic regression is introduced to learn models for the nonlinear stressstrain relation. Our aim is to identify models given a set of different input features without restricting the search by too many modelling assumptions.

II. Nonlinear constitutive stress-strain relation

The challenge in RANS-based turbulence modelling is the closure of the RANS equations by a model for the Reynolds-stress tensor τ i j . This symmetric second order tensor can be decomposed into

τ i j = 2k (a i j + 1 3 δ i j ), (1)
in which a i j represents the non-dimensionalised anisotropic part and k the turbulent kinetic energy the isotropic part.

Based on the decomposition of the mean velocity gradient tensor into the symmetric and the anti-symmetric part

∂ j U i = s i j + ω i j , (2)
with the mean strain-rate tensor s i j = 1 2 (∂ j U i + ∂ i U j ) and the rotation rate tensor

ω i j = 1 2 (∂ j U i -∂ i U j ),
Pope was the first to use the Cayley-Hamilton theorem in order to derive an integrity basis of ten nonlinear base tensors and five corresponding invariants [START_REF] Pope | A more general effective-viscosity hypothesis[END_REF]. Given a turbulent time scale τ to non-dimensionalise S i j = τs i j and Ω i j = τω i j the base tensors can be combined linearly as a functional form for a i j :

a i j (S i j , Ω i j ) = N n=1 T (n) i j α n . (3)
The base tensors are

T 1 i j = S i j , T 2 i j = S ik Ω k j = Ω ik S k j , T 3 i j = S ik S k j - 1 3 δ i j S mn S nm , T 4 i j = Ω ik Ω k j - 1 3 δ i j Ω mn Ω nm , T 5 i j = Ω ik S k l S l j -S ik S kl Ω l j , T 6 i j = Ω ik Ω k l S l j + S ik Ω kl Ω l j - 2 3 S kl Ω l m Ω mk , T 7 i j = Ω ik S k m Ω mn Ω n j -Ω ik Ω k m S mn Ω n j , T 8 i j = S ik Ω k m S mn S n j -S ik S k m Ω mn S n j , T 9 i j = Ω ik Ω k m S mn S n j + S ik S k m Ω mn Ω n j - 2 3 S k m S mn Ω n p Ω pk , T 10 i j = Ω ik S k m S mn Ω n p Ω p j -Ω ik Ω k m S mn S n p Ω p j (4)
and the corresponding invariants read

I 1 = S mn S nm , I 2 = Ω mn Ω nm , I 3 = S k m S mn S nk I 4 = Ω k m Ω mn S nk , I 5 = Ω k m Ω mn S n p S pk . (5)
While the combination of base tensors is linear, according to Pope's analysis the scalar coefficient fields α n are nonlinear functions of the invariants, i.e. α n = α n (I 1 , ..., I 5 ). The identification of the functional form of the coefficients is the essential step to build a nonlinear eddy-viscosity model. Classical methods to identify the functional forms are given in [START_REF] Pope | A more general effective-viscosity hypothesis[END_REF][START_REF] Apsley | A New Low-Re Non-Linear Two-Equation Turbulence Model for Complex Flows[END_REF]. As mentioned above these methods are based on modelling assumptions such as the weak equilibrium hypothesis, which omits the transport of anisotropy [START_REF] Gatski | Nonlinear eddy viscosity and algebraic stress models for solving complex turbulent flows[END_REF]. In the following, optimal scalar fields α n are identified for the given test cases and further used to learn a functional form for α n using symbolic regression.

III. Test case, high-fidelity data and solver

The test case, for which we want to infer and test correction models, is the flow over a series of periodic hills in 2D [START_REF] Mellen | Large Eddy Simulation of the flow over periodic hills[END_REF], the geometry is shown in Fig. 1. The flow exhibits separation on the leeward side of the hills and reattachment in between. Both features are challenging to predict by LEV models [START_REF] Jakirlic | Extended excerpt related to the test case: "Flow over a periodical arrangement of 2D hills[END_REF]. We consider two different Reynolds-numbers Re H = 2800 and 10595 based on the bulk-velocity at the hill's crest and the hill height. At both the lower as well as the upper wall a no-slip boundary condition is applied. The flow setup is periodic in the stream-wise direction and driven with a volume forcing. The CFD solver used is simpleFoam from the OpenFOAM toolbox [START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF]. The inferred correction models resulting from the proposed learning method are implemented in a new turbulence model based on the k -ω model. The high-fidelity data used for model-learning and validation are mean velocity fields and Reynolds-stress fields from Breuer et al. [START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Reynolds numbers[END_REF]. More precisely, DNS data for the lower Reynolds number case and LES data for the higher Reynolds number one.

IV. Optimal coefficients for the nonlinear constitutive relation

As the modelling of the anisotropy of the Reynolds-stress tensor using the linear stress-strain relation is insufficient, we introduce an additive term b ∆ i j to account for the model-form error

τ i j = -2ν t S i j + 2 3 kδ i j + 2kb ∆ i j , (6)
which is a non-dimensional symmetric tensor with zero trace. The mean strain-rate S i j , the Reynolds-stress tensor τ i j and the turbulent kinetic energy k are directly available from high-fidelity simulations (LES or DNS) or from experiments. The eddy viscosity ν t can be obtained by passively solving the turbulent transport equations of a turbulence model, e.g. k -ω, for k and the mean velocity field U i (frozen approach) [START_REF] Weatheritt | The development of algebraic stress models using a novel evolutionary algorithm[END_REF]. The model-form error b ∆ i j is then readily computable as a residual given the equation above.

The obtained model-form error is case-dependent, meaning that it only corrects the stress-strain relation for a given setup (geometry, boundary conditions, Reynolds number, etc.). The key is to distill information from it to identify a suitable mathematical expression as a corrective term, which can be used for prediction of other flow cases. For that we identify the coefficient fields of the base tensor series given input data b ∆ i j . The relation between the data and the base tensor series is expressed by a statistical model

b ∆ = H a + , (7)
in which the sparse block-diagonal matrix H ∈ R 6K ×M with M = N × K contains the base tensors ordered per component and spatial location k = 1, ..., K

H =                                       T (1) x x,1 • • • T (N ) x x,1 0 • • • 0 . . . 0 • • • 0 T (1) x y,1 • • • T (N ) x y,1 0 • • • 0 . . . 0 • • • 0 T (1)
x z,1

• • • T (N ) x z,1 0 • • • 0 . . . 0 • • • 0 T (1) y y,1 • • • T (N ) y y,1 0 • • • 0 . . . 0 • • • 0 T (1) y z,1 • • • T (N ) y z,1 0 • • • 0 . . . 0 • • • 0 T (1) z z,1 • • • T (N ) z z,1 0 • • • 0 . . . 0 • • • 0 0 • • • 0 T (1) x x,2 • • • T (N ) x x,2 . . . 0 • • • 0 . . . . . . . . . . . . . . . . . . 0 • • • 0 T (1) z z,2 • • • T (N ) z z,2 . . . 0 • • • 0 . . . . . . . . . . . . . . . . . . . . . 0 . . . 0 0 . . . 0 . . . T (1) z z, K • • • T (N ) z z, K                                       (8)
and the vectors a ∈ R M and b ∆ ∈ R 6K contain the stacked coefficient and tensor fields as well ordered per component and spatial location, respectively:

a = (α 1,1 , . . . , α N,1 , α 1,2 , . . . , α N, K ) T , b ∆ = (b ∆ x x,1 , b ∆ x y,1 , b ∆ xz,1 , b ∆ y y,1 , b ∆ y z,1 , b ∆ z z,1 , b ∆ x x,2 , . . . , b ∆ z z, K ) T . (9)
The random variable serves as an additive error term expressing the fact, that given the best set of coefficients α n,k the base tensor series might not fit the data perfectly. As the input data (Reynolds-stress and mean velocity field) is relatively smooth spatially, we also want to infer smooth coefficients fields. In a preliminary study utilising ordinary least-square regression, we discovered rough coefficient fields, in which the values next to each other differed several orders of magnitude, which is a common issue for inverse problems [START_REF] Hansen | Rank-Deficient and Discrete Ill-Posed Problems[END_REF]. While these coefficients lead to a very low reconstruction error of the anisotropic Reynolds-stress, they are discarded as they are not useful to serve as a target for the symbolic regression. A lack of smoothness leads to unreasonably complex models, with an increased risk of overfitting the rough coefficient fields. Furthermore, more complex correction models also increase the numerical stiffness of the CFD solver leading to a reduction or loss of numerical stability. Even though the smoothness constraint increases the a priori reconstruction error, it enables the discovery of sparse correction models.

The inference of the coefficients given the data b ∆ is done via Bayesian inversion [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF][START_REF] Tarantola | Inverse Problem Theory[END_REF] p(a| b ∆ ) ∝ L(b ∆ |a) p(a) [START_REF] Mcconaghy | FFX: Fast, Scalable, Deterministic Symbolic Regression Technology[END_REF] in which the posterior probability p(a|b ∆ ) is the product of the likelihood L(b ∆ | a) based on the statistical model defined above and the prior probability p(a), which embodies the smoothness constraint. As the dimension of this problem depends on the mesh size, which is already large for the considered 2D flow case (18000 cells), we focus on the inference of the maximum a priori estimation (MAP) of the posterior

a M AP = arg max a p(a| b ∆ ) (11)
and not the full posterior distribution. We further simplify the problem by assuming that both the likelihood function and the prior are normally distributed

L(b ∆ |a) ∝ exp - 1 2 (b ∆ -H a) T R -1 (b ∆ -H a) , (12)
p(a) ∝ exp - 1 2 (a -a o ) T P -1 (a -a o ) . (13)
Thus, also the posterior is normally distributed and, under the assumption of a prior zero mean a o = 0, it reads

p(a| b ∆ ) ∝ exp - 1 2 (a -a M AP ) T Σ -1 (a -a M AP ) (14)
with mean a M AP and covariance matrix Σ

a M AP = K b ∆ , ( 15
) Σ = (I -K H )P, (16)
K = PH T (R + H PH T ) -1 ( 17
)
the latter is known as the Kalman gain [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. In the present study the observation covariance matrix R ∈ R 6K is assumed to be diagonal R = σ R I implying that the underlying spatial correlation of the data b ∆ is omitted. Further work will also use other covariance structures. The prior covariance P ∈ R M acts to regularise the inversion problem in order to obey the smoothness constraint on the coefficient fields and is defined as

P(x, x ) =        σ P exp -(x-x ) 2 L 2 , if exp -(x-x ) 2 L 2 ≥ 0.001 0, otherwise, (18)
in which σ P is the scalar variance, x represents the coordinates of location within the mesh and L is the correlation length. The Gaussian kernel is a common choice resulting in a smooth correlation function. With a correlation length shorter then the smallest cell size the MAP-inversion framework reduces to a regularised least-squares problem, while a larger correlations length increases both the smoothness of the coefficients and the reconstruction error.

Given values for σ R , σ P and L we are now able to compute a M AP with [START_REF] Gatski | Nonlinear eddy viscosity and algebraic stress models for solving complex turbulent flows[END_REF]. As the size of the problem and therefore the size of the matrices present in the Kalman gain in [START_REF] Jakirlic | Extended excerpt related to the test case: "Flow over a periodical arrangement of 2D hills[END_REF] scales with the number of mesh points of the flow case at hand, we need to exploit the sparsity of the matrices in order to store them efficiently and make the computations tractable. Due to their structure defined above especially the memory requirements for H and R can be drastically reduced. The density of P given [START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF] is controlled by a cut-off for small-sized correlations. As it depends on the correlation length L, also its memory requirements increase with a larger correlations length. For the periodic hill test case we have tested L = 0.01 and L = 0.001, see Fig. 2a and 2b respectively. In these figures, the reconstruction error

= (b ∆ -H â) T (b ∆ -H â) 0.5
for an inferred set of coefficients â is shown against an empirical roughness criterion

s = (|∇ â0 |) max (|∇ â0 |) mean , (19)
corresponding to the maximum of the gradient of â0 = α0,k normalised by its mean value. Unregularised inversion is prone to small differences in the data leading to large discrepancies in the inferred coefficients [START_REF] Hansen | Rank-Deficient and Discrete Ill-Posed Problems[END_REF] , i.e. localised large gradients. The rationale of the criterion is therefore to identify if unphysically large gradients are present for a given regularisation strategy. In general, the lower s is, the lower are the highest gradients, which contributes to smoother coefficient fields. The coloured marker-lines represent a different number of employed base tensors ranging from N = 2 to 5, which is a common choice for EARSM for 2D [START_REF] Leschziner | Statistical Turbulence Modelling for Fluid Dynamics -Demystified: An Introductory Text for Graduate Engineering Students[END_REF], and the scalar variances of P and R are summarised to

λ = σ R σ P , (20)
equivalent to a Tikhonov regularisation parameter in case of regularised least-square regression [START_REF] Hansen | Rank-Deficient and Discrete Ill-Posed Problems[END_REF]. In general, the lower the regularisation parameter λ, the lower the reconstruction error and the larger s, and vice versa. For λ ≤ 0.001 the reconstruction error decreases while the roughness increases more than one order of magnitude. By increasing the number of base tensors from N = 2 to N = 3 the reduction of the reconstruction error is significant for a given λ, but only minor for a further expansion. This implies that for a correction model at least N = 3 should be considered. For λ ≥ 0.001 the reconstruction error increases but the smoothness settles, also the number of base tensors is not important anymore in terms of the reconstruction error, meaning that the problem is dominated by the regularisation.

Changing the correlation length (L = 0.01 in Fig. 2a, L = 0.001 in Fig. 2b) decreases the roughness of the coefficient fields slightly, but has a minor effect overall. We have observed for even larger correlation lengths a further decrease in roughness for a given λ, but also an increase in the reconstruction error. Smoothing the coefficient fields by increasing L is therefore detrimental to the endeavour of designing a prior to steer the inversion towards low roughness and low reconstruction error (lower-left corner of the Fig. 2).

Finally, the preceding parametric study helps us in identifying suitable coefficient fields to serve as targets for the symbolic regression. We identify coefficient fields given L = 0.01 and λ = 0.0001 corresponding to a reasonable compromise between low bias according to a low reconstruction error and low roughness s. Further research will focus on a rigorous optimisation of the parameters λ and L to find an optimum of low reconstruction error and low roughness.

V. Deterministic Symbolic Regression of Correction Models

This Section deals with the methodology of deterministic symbolic regression. First, it is explained how a library of candidate functions is build automatically given raw input features and a set of mathematical operations. Second, details are given on the algorithm to identify the best linear combination of as few as possible candidate functions to fit the target coefficient fields, which were identified in the preceding section.

A. Building a library of candidate functions

The deterministic symbolic regression constructs an over-complete library of possible nonlinear candidate functions B approximating given data and identifies the relevant ones by adopting a sparsity constrain. An important step is the design of the library of candidate functions. Because the correction models we want to identify are based on the nonlinear eddy viscosity concept, we use the invariants I 1 and I 2 of the nonlinear base tensor expansion as our primitive input features, which are shown for the periodic hill test case calculated from the mean velocity of the LES data of [START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Reynolds numbers[END_REF] at Re H = 10, 595 in Fig. 3a and3b respectively. Further candidate functions are constructed by applying mathematical operations to the already existing ones present in the library of candidate functions, see Table 1. The operations can be applied in a different order and also repetitively. For example a library containing only polynomials and products of the resulting set of candidates is encoded in

B = {R | P | M} , (21)
in which the vertical line | indicates that the candidate functions resulting from the operation to the left are passed on to the operation to the right. The resulting new functions after each operation are appended to the library. This procedure makes the construction of the library automatic similar to FFX [START_REF] Mcconaghy | FFX: Fast, Scalable, Deterministic Symbolic Regression Technology[END_REF]. 1 Operations to build a library of candidate functions Since the procedure can produce duplicates at different stages of the process, we check at the end for equivalent expressions and retrieve only a unique set of candidate functions. In addition for every candidate function it is checked, if it can be evaluated for the flow case, for which we want to find a correction model. Functions, which are not defined for the input domain provided by the data, are deleted from the library. Obviously the nonlinearity of the resulting library B can be determined by the type and order of the operations and by their frequency. This is illustrated in Table 2.

Symbol Operation Details

R (•) raw input features P (•) p p = [±0.5, ±2, ...] F f (•) f = [sin, cos, exp, log, ...] A (•) + (•) summation M (•) × (•) product Table

Order

example of resulting most complex candidates 2 Illustration of resulting most complex candidate functions for a given order of operations.

R | P | M I 0.5 1 I 2 2 R | P | F | M exp(I 0.5 1 ) log(I 2 2 ) R | P | A | F exp(I 0.5 1 + I 2 2 ) Table
As a starting point for the present work we focus on a weakly nonlinear library using the rule described in [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF] with p = {1.0, 2.0}. In addition to the two invariants I 1 and I 2 we also include a constant function c to the set of raw input features. This helps to identify very sparse models, for which the form of a base tensors only needs to be scaled by a factor and not altered by a spatially-dependent function. The resulting library B has 16 candidate functions and reads

B = c, I 1 , I 2 , I 2 1 , I 2 2 , I 2 1 I 3 2 , I 4 1 I 2 2 , I 1 I 2 2 , I 1 I 3 2 , I 1 I 4 2 , I 3 1 I 2 , I 2 1 I 4 2 , I 2 1 I 2 , I 1 I 2 , I 3 1 I 2 2 , I 2 1 I 2 2 (22)
which are evaluated using the reference data. Thus, by using p = {1.0, 2.0} and allowing products of candidate function with the same base, we effectively achieve exponents p = {1.0, 2.0, 3.0, 4.0}. The evaluated candidate functions are stored column-wise in the library matrix

B =          c| k=0 I 1 | k=0 . . . I 2 1 I 2 2 | k=0 . . . . . . . . . c| k=K I 1 | k=K . . . I 2 1 I 2 2 | k=K          . ( 23
)

B. Deterministic Symbolic Regression

The deterministic symbolic regression identifies the most relevant candidate functions by solving the optimisation problem

Θ (n) = arg min Θ(n) B Θ(n) -a n 2 2 + λ r Θ(n) q , (24)
in which the vector of coefficients Θ (n) needs to be identified. The target is a specific coefficient field a n for a given base tensor n. The regularisation term using norm q = 1 (lasso regression) or q = 2 (ridge regression) acts to increase the sparsity of Θ (n) , i.e. increasing the number of zeros in order to turn off the corresponding candidate functions [START_REF] Brunton | Discovering governing equations from data: Sparse identification of nonlinear dynamical systems[END_REF][START_REF] Rudy | Data-driven discovery of partial differential equations[END_REF], proportional to the magnitude of the Tikhonov parameter λ r . We use STRidge (Sequential Threshold Ridge Regression), which solves the optimisation problem in Eq. 24 using ridge regression q = 2 iteratively [START_REF] Brunton | Discovering governing equations from data: Sparse identification of nonlinear dynamical systems[END_REF]. After each step sparsity is achieved by setting coefficients to zero, which have a smaller magnitude than a given threshold t max , and the regression is repeated for the non-zero coefficients until a suitable model is identified. For the resulting set of coefficients an ordinary least-square regression, i.e. λ r = 0, is performed to achieve the unbiased values for the coefficients. Once Θ (n) has been identified for every base tensor n, the mathematical expression of the discovered correction model can be retrieved from

M := b ∆ i j = N n=1 BΘ (n) T (n) i j . (25)
Using data the flow over periodic hills at Re H = 10595 [START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Reynolds numbers[END_REF] only one step of STRidge with λ r = 0.02 was performed and the threshold was defined as t max = ξ max(| θ|) with 0 < ξ < 1 depending on the largest absolute coefficient max(| θ|). This variation of STRidge allows to determine a hierarchical set of models given different values of ξ. Give ξ = 0.9 the expression for the resulting correction model reads M (1) = 7.0815 I 1 T (1) i j + 4.2099 T (2) i j + 3.6909 I 1 T (3) i j -0.0509 T (4) i j

(26)

and by applying the same setup for the symbolic regression to the data of Re H = 2800 we obtain M (2) = -6.4069 I 2 T (1) i j + 4.1292 T (2) i j -17.4622 I 2 T (3) i j + 18.1459 I 2 T (4) i j .

The models have a very sparse structure as only the raw input features and the constant have been identified from B to serve as functions for the coefficients. This is a result of a large ξ = 0.9 value. Interestingly, for M (1) only I 1 with a positive coefficient and for M (2) only I 2 with a negative coefficient have been chosen for T (1) i j and T (3) i j . This choice is reasonable, as the invariants exhibit a similar spatial structure with opposite signs shown in Fig. 3. While the coefficient for T (3) i j is almost the same for both models, the identified function for T (4) i j is different. The a priori error is calculated using two metrics, the root-mean squared error and the tensor alignment ρ between the training data and the model

(M) = 1 6K K k=1 3 i=1 i j=1 M i j, k -b ∆ i j, k 2 (28) ρ(M) = K k=1 M i j, k b ∆ i j, k M mn, k M nm, k b ∆ pq, k b ∆ q p, k , (29)
in which the tensors components i j are evaluated at each grid point k. The results are given in Table 3. Given the l2-norm of the distance between the model and the data in eq. 24 the metric quantifies the achieved goodness-of-fit of the symbolic regression. In addition, the metric ρ calculates how well the resulting tensor shape is aligned with the data, and therefore does not include magnitude information. In comparison to the literature the achieved a priori error is similar to other approaches using deep learning [START_REF] Ling | Reynolds averaged turbulence modelling using deep neural networks with embedded invariance[END_REF] or genetic programming based symbolic regression [START_REF] Weatheritt | A novel evolutionary algorithm applied to algebraic modifications of the RANS stress-strain relationship[END_REF].

Model Training case λ L λ r ξ ρ M (1) Re H = 10595 0.0001 0.01 0.02 0.9 0.0875 0.8197 M (2) Re H = 2800 0.0001 0.01 0.02 0.9 0.0768 0.7981 Table 3 Parameters in derivation and a priori error of the models.

VI. Prediction

Fig. 4 shows predictions of the stream-wise velocity for Reynolds number Re H = 10595 using the correction models M (1) and M (2) . While M (1) has been derived from the same dataset, using it to predict the same flow serves as a verification. In general, both models show an improvement compared to the baseline LEVM k -ω all over the domain. The most challenging quantity of this flow case is the reattachment point x a between the hills. This is reflected by a vanishing wall-shear stress at the lower wall, shown in Fig. 5 and summarised in Tab. 4. Again the predictions of the models are closer to the reference data than the baseline LEVM, which expresses the fact that the symbolic regression successfully identified corrections of the baseline LEVM. Also for all shown quantities M (1) is closer to the reference data than M (2) consistently. However, for the lower Reynolds-number case (Re H = 2800) shown in Fig. 6, M (1) performs similarly well and predicts the reattachment closer to the reference than M (2) . This behaviour is not unexpected as the models are derived to be sparse corrections of the underlying LEVM in order to avoid over-fitting. Further work needs to be done in order to study the dependency of the sparsity-controlling cut-off parameter ξ on the predictive performance of the resulting models in relation to the Reynolds-number. Fig. 4 Predictions of the stream-wise velocity for Re H = 10595 using correction model M (1) , k -ω LEVM M o and EARSM M E ARS M from literature [START_REF] Hellsten | New Advanced k-w Turbulence Model for High-Lift Aerodynamics[END_REF]. Validation data from [START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Reynolds numbers[END_REF].

For the present case, the Explicit-Algebraic Reynolds-stress model (EARSM) from literature [START_REF] Hellsten | New Advanced k-w Turbulence Model for High-Lift Aerodynamics[END_REF] performs worse than the linear k -ω model all over the domain at Re H = 10595, see Fig. 4. In a benchmark study [START_REF] Jakirlic | Extended excerpt related to the test case: "Flow over a periodical arrangement of 2D hills[END_REF], in which several EARSM and RSM were applied to this test case, it was also observed that not necessarily more complex models (1) , k -ω LEVM M o and EARSM M E ARS M from literature [START_REF] Hellsten | New Advanced k-w Turbulence Model for High-Lift Aerodynamics[END_REF]. Validation data from [START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Reynolds numbers[END_REF]. perform better.

Re H

M (1) M (2) M o LES [START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Reynolds numbers[END_REF] x a 2800 5.69 6.17 7.61 5.37 x a 10595 5.02 5.42 5.96 4.54 Table 4 Predictions of reattachment point x a for Re H = 10595 using correction model M (1) , M (2) and k -ω LEVM M o in comparison to the LES data from [START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Reynolds numbers[END_REF].

VII. Conclusion

In this contribution we show that symbolic regression is a powerful machine learning method to learn the mathematical form of models from data for the purpose of RANS-based turbulence modelling. The goal is to identify concise models ready to be implemented in existing codes as used in industry proven in [START_REF] Weatheritt | A novel evolutionary algorithm applied to algebraic modifications of the RANS stress-strain relationship[END_REF][START_REF] Weatheritt | The development of algebraic stress models using a novel evolutionary algorithm[END_REF] but with a deterministic framework, which can be applied to high-dimensional regression problems.

We introduced a framework for the inference of optimal coefficients fields of the base tensor series used in the nonlinear eddy viscosity concept. By utilising a maximum a posteriori (MAP) based Bayesian inversion we conduced a parameter study in order to achieve coefficients with low roughness and low reconstruction error by using data from high-fidelity simulations of the flow over periodic hills at Re H = 2800 and 10595.

Fig. 1

 1 Fig. 1 Stream-wise velocity of flow over periodic hills in 2D at Re H = 10595. LES data from [19].

Fig. 2

 2 Fig. 2 Reconstruction error versus roughness metric s for two different correlation length L. The number of base tensors N are indicated by coloured marker lines. The regularisation parameter λ are associated to the points next to the mentioned value.

2 Fig. 3

 23 Fig.3Two invariants of the nonlinear base tensor series for flow over periodics hills at Re = 10, 595. Data from[START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Reynolds numbers[END_REF].

Fig. 5

 5 Fig.5Prediction of the wall-shear stress for Re H = 10595 using correction model M(1) , k -ω LEVM M o and EARSM M E ARS M from literature[START_REF] Hellsten | New Advanced k-w Turbulence Model for High-Lift Aerodynamics[END_REF]. Validation data from[START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Reynolds numbers[END_REF].

Fig. 6 Predictions of the stream-wise velocity for Reynolds number Re H = 2800 using two correction models M (1) and M (2) and k -ω LEVM M o . Validation data from [START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Reynolds numbers[END_REF].

The resulting coefficient fields were used as targets of the symbolic regression based on sparsity-promoting regularised least square regression. We have successfully identified models and presented the a priori error and the predictive performance of both of them.

The method shows potential for data-driven correction of RANS-based turbulence modelling. Next steps will be the application to flow cases with a different geometry and establishing a relation between the a priori error and the predictive performance of a model. In addition, the two step procedure of inferring the targets and conducting the symbolic regression separately will be merged into a single optimisation framework.