Direct Numerical Simulation of hydrogen combustion at auto-ignitive conditions: Ignition, stability and turbulent reaction-front velocity

Direct Numerical Simulations (DNS) are performed to investigate the process of spontaneous ignition of hydrogen flames at laminar, turbulent, adiabatic and non-adiabatic conditions. Mixtures of hydrogen and vitiated air at temperatures representing gas-turbine reheat combustion are considered. Adiabatic spontaneous ignition processes are investigated first, providing a quantitative characterization of stable and unstable flames. Results indicate that, in hydrogen reheat combustion, compressibility effects play a key role in flame stability and that unstable ignition and combustion are consistently encountered for reactant temperatures close to the mixture’s characteristic crossover temperature. Furthermore, it is also found that the characterization of the adiabatic processes is also valid in the presence of non-adiabaticity due to wall heat-loss. Finally, a quantitative characterization of the instantaneous fuel consumption rate within the reaction front is obtained and of its ability, at auto-ignitive conditions, to advance against the approaching turbulent flow of the reactants, for a range of different turbulence intensities, temperatures and pressure levels. © 2021 The Author(s). Published by Elsevier Inc. on behalf of The Combustion Institute. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Hydrogen-firing of stationary gas turbines is emerging as one of the most robust approaches to reduce carbon emissions from large-scale power generation. This equally applies, in a convenient synergy, to power generation schemes that can utilize a steady stream of hydrogen from large-scale reforming of natural gas with carbon capture and storage (CCS) [1] or, exploiting excess power from non-dispatchable renewable energy resources (wind and solar), an unsteady stream of hydrogen produced from water electrolysis coupled to large-scale energy storage solutions (power-to-H 2 -to-power) [2] .
However, state-of-the-art gas-turbine technology does not presently allow, without important performance compromises, for combustion of pure (undiluted) hydrogen. This fuel notoriously poses important burner design challenges with respect to flame stability and NO x emissions that are conventionally solved by dilution of the hydrogen fuel with large quantities of steam or nitrogen [3] . The main reason for these problems is due to hydrogen's higher reactivity compared to natural gas, the standard gaseous fuel for gas turbines. Hydrogen's high reactivity introduces severe challenges in simultaneously achieving low emission performance together with static and dynamic flame stability (i.e. avoiding flashback [4] and thermo-acoustics instabilities [5,6] ), and remains one of the main obstacles for large-scale, clean and efficient utilization of hydrogen in gas turbines.
In this context, multi-stage combustion systems seem to offer the most promising solution for power plants that, in today's changing power market, have to ensure high turndown ratios, part-load efficiency, and fuel flexibility (including hydrogen firing), while keeping pollutant emissions low. A two-stage premixed system, in which the two combustion stages are distributed longitudinally in a sequential arrangement and separated by an air-dilution section is presently employed in the Ansaldo GT36 Hclass gas turbine [7] . Here, the first combustion stage consists of a conventional aerodynamically/propagation-stabilized flame and, if operated alone, provides optimal part load performance (high efficiency, low emission) while the second combustion stage consists of a so-called reheat flame stabilized (mainly) by spontaneous ignition in a sequential combustor [8] . In Ansaldo's longitudinallystaged combustion system, the first combustion stage serves as a hot-gas generator while the predominant energy conversion occurs in the second stage, posing new and interesting challenges due to the unconventional combustion conditions and rate-limiting processes that characterize these reactive flows. Other gas turbine manufacturers are pursuing similar longitudinal fuel staging strategies although these are typically characterized by different staging ratios, in which most of the fuel is consumed in the first stage, ultimately resulting in different combustion and operational behavior [9,10] .
In principle, the reheat combustion scheme, due to its reliance on spontaneous ignition to achieve flame stabilization [11] , is well-suited to provide intrinsically stable and clean combustion of hydrogen-rich fuel mixtures, as recently demonstrated [12] . This result is achieved through an operational strategy for hydrogenfiring that implements a reduction of the flame temperature in the first stage, through lean(er) operation of the propagation-stabilized flame, achieving flashback and NO x control locally and, simultaneously, ensuring a beneficial increase in the ignition delay time of the (now colder) reactants' mixture entering the second stage. There, in the sequential combustor, the reactants' inlet temperature represents the main rate-controlling parameter that controls, to leading order, the stabilization location of the reheat flame through the process of spontaneous ignition. Therefore, this is one of the key quantities focused upon in the present paper whereas the fueloxidant equivalence ratio plays itself a minor role with respect to flame stability (within the relevant operational range of interest). The reliance on spontaneous ignition rather than flame propagation to achieve stabilization of the flame in the second (main) combustion stage has several advantages. The principal consequence of this strategy is that high bulk velocities can be utilized within the sequential combustor flow path to reduce NO x formation (due to short residence time) and diminish the propensity for flashback (due to high flow velocity). Moreover, increasing the fuel supply in the second stage, while not affecting flame stabilization (principally controlled by the reactant's temperature), fully compensates for the reduced fuel addition in the first stage, maintaining the target flame temperature in the second stage and minimizing or eliminating de-rating of the engine [12] .
Depending on the boundary and operating conditions in practical combustion applications, it is reasonable to expect that complex mixed combustion modes can occur in the sequential combustor [13] . These are characterized by the simultaneous existence of deflagration and spontaneous ignition fronts, either present in different localities or co-located, and their interaction makes the physical process of reheat combustion more complex to understand and predict. The local behaviour of propagating flames or spontaneous auto-igniting fronts is affected by their surroundings. As a result, the balance of combustion modes is able to affect global combustor behaviour through feedback mechanisms with the velocity and acoustic fields (leading to flashback or to thermoacoustic instabilities).
There have been relatively few past studies on combustion under reheat conditions, i.e. vitiated oxidant, elevated pressures (up to ∼ 25 bar ) and high reactant temperatures ( > 10 0 0 K ). One of the early research efforts was conducted by the Institute of Combustion Research of DLR. Pressurized laboratory experiments were performed on a scaled, geometrically simplified version of the sequential combustor fired with hydrocarbon fuels [14][15][16] . How-ever, these first experimental studies, focused on the formation of auto-ignition kernels in the mixing duct, offer limited insight about the reheat combustion process in the main flame stabilization region. More recently, joint numerical and experimental studies were performed at ETH and PSI (Zurich, Switzerland). The researchers performed state-of-the-art Large Eddy Simulation (LES) with the Dynamically Thickened Flame [17] to investigate the occurrence of deflagration and spontaneous ignition in methane-air flames and their relative importance in flame stabilization [18][19][20][21][22] . In the numerical modelling study by Krisman [23] , the designation of a uniquely-defined, quantitative reference speed for laminar premixed flames at auto-ignitive condition has been proposed for the first time. Other modelling studies specifically focused on the dynamic response of auto-igniting flames, using LES to extract the flame transfer function (FTF) based on pressure and velocity fluctuations [24,25] while the most recent investigations have highlighted the importance of inlet temperature fluctuations that must be accounted for in a 3 ×3 flame transfer matrix (FTM) [26][27][28] .
Beyond the aforementioned efforts mostly focusing on hydrocarbon fuels, to date, only a handful of studies have investigated the characteristic features of hydrogen combustion at reheat conditions. The earliest among these were zero-dimensional (0D) and one-dimensional (1D) reactor modelling studies characterizing ignition and propagation time scales, complementing fullscale, high-pressure experiments [29][30][31] . Only very recently fullfledged, three-dimensional Direct Numerical Simulations (DNS) of turbulent premixed hydrogen-air combustion at reheat conditions (albeit atmospheric pressure) have been performed in conjunction with detailed chemical kinetics and Chemical Explosive Mode Analysis [32] to quantify the relative importance of flame propagation versus spontaneous ignition for a range of turbulence intensities in statistically planar flames [33] and in the presence of wall heat loss in a semi-realistic combustor geometry [34] .
The present effort builds upon the aforementioned DNS studies and deploys one-two-and three-dimensional DNS to investigate the conditions leading to steady or unsteady ignition and combustion in premixed hydrogen-air reheat flames under laminar, turbulent, adiabatic and non-adiabatic conditions. The first part of this study focuses on the unsteadiness related to the spontaneous ignition process itself that takes place when combustion arises in a mixture of preheated reactants. The occurrence of selfexcited flame instabilities, emerging well after their initial ignition and induced by a variation of the reactant temperature, is reported ( Section 3 ). The second part of the present study reports the effect of the turbulence intensity characterizing the approaching reactant flow on the turbulent reheat flame-brush mean displacement velocity ( Section 4 ). Section 2 briefly describes the physical process that is the objective of this investigation and the numerical tool chosen for the study (Sandia's S3D DNS code) while the main results are summarized in Section 5 , which also presents an outlook about further work on the topic.

Spontaneous ignition of hydrogen reheat flames
For hydrogen-rich reactant mixtures, the ignition delay (or induction time), τ ig , is a key quantity at preheated reheat conditions that is strongly dependent on initial temperature ( T 0 ), pressure ( P 0 ) and composition ( Y k, 0 ). The ignition delay exhibits a large rate of change near the crossover temperature where the reaction rate of the elementary H + O 2 branching step equals that of the recombination step H + O 2 + M [35] . This characteristic behavior of hydrogen-air systems has been extensively investigated in the past and it can be accurately reproduced through advanced and elegantly simple kinetic models [36] . It has been observed that, when transitioning from methane-air mixtures to hydrogenenriched methane-air mixtures and finally to hydrogen-air mixtures, the temperature dependence of τ ig at conditions relevant to reheat combustion increases considerably with hydrogen content [37] . Therefore, for hydrogen-air mixtures, in which the fuel consists of pure hydrogen undiluted by hydrocarbons, nitrogen or steam, it is reasonable to expect a spontaneous ignition behavior that largely departs from the behavior typically observed in hydrocarbon-air premixed flames.
Nearly 40 years ago Zel'dovich postulated the existence of a "spontaneous propagation" regime of relevance for combustion at auto-ignitive conditions that is characterized by an intermediate ignition-front velocity between the deflagration and the detonation velocities [38] . Spontaneous propagation of the ignition front occurs if the inverse of the magnitude of the local ignition-time gradient is larger than the deflagration velocity and smaller than the detonation velocity, i.e.
where c is the speed of sound and S L and S D are the "conventional" deflagration and detonation velocities, respectively. Based on Zel'dovich criterion, the transition between deflagration and spontaneous propagation of ignition fronts is highly sensitive to spatial gradients in temperature present in the reactant mixture and to the temperature dependence of the ignition delay. Consequently, it is reasonable to assume that compressibility effects (e.g. compression heating) can play a key role in controlling the behavior of reactive flows in the spontaneous-propagation regime. The above assumption has been already validated numerically through previous DNS studies, focusing on compressionignition engine conditions [39,40] , that illustrated a dependency between the propagation speed of the front and thermal-or combined thermal-and composition-gradients, modulated by turbulent mixing and isentropic compression heating. These early observations are expected to be of particular relevance to hydrogen-air systems due to the large values that characterize d τ ig /d T 0 in temperature ranges near crossover.

Direct Numerical Simulation code
From the aforementioned discussion it is important to capture the coupling between pressure, density and temperature fields for the reactive flows under investigation, in addition to a detailed representation of the chemical reaction kinetics. To this end, the compressible reacting Direct Numerical Simulation (DNS) code, S3D, originally developed at Sandia National Laboratories [41] , is employed for all calculations described in Sections 3 and 4 .
S3D is written in FORTRAN 90 and uses the Message Passing Interface (MPI) for interprocess communication in parallel execution. In the present application, the algorithm implemented in S3D solves the Navier-Stokes equations for a compressible fluid in conservative form on structured, Cartesian meshes in one-, two-and three-dimensional computational domains to simulate premixed combustion of H 2 -air flames at auto-ignitive conditions representative of a reheat combustion system. The present simulations are mostly limited to atmospheric pressure due to computational cost with increasing pressure. All DNS use a spatial resolution of at least δs = δx = δy = δz = 25 μm ( δs = 10 μm in 1-D and 2-D calculations) that is sufficient to resolve all spatial scales of the reactive flows investigated at atmospheric pressure. The spatial derivatives are computed with an eighth-order, explicit, centered, finite-difference scheme (third-order one-sided stencils are used at the domain boundaries in the non-homogeneous directions) in conjunction with a tenth-order, explicit, spatial filter, as suggested by Kennedy and Carpenter [42] , to remove high frequency noise and reduce aliasing error. A fourth-order, six-stage, explicit Runge-Kutta scheme, described in [43] , is used for time integration and the time step is set to δt = 4 ns for all reactive flows investigated.
Thermodynamic properties are modelled as polynomial functions of temperature and transport coefficients as described in the CHEMKIN and TRANSPORT packages, respectively [44] . Radiative heat transfer is neglected because of the modest optical thickness of hydrogen flames. The chemical reactions in the gas phase are described by a detailed mechanism for hydrogen combustion in air [45] . This mechanism consists of 9 species and 19 elementary reaction steps. Nitrogen is assumed to be inert such that NO xformation reactions are not considered.
Inflow and outflow boundary conditions are implemented following the Navier-Stokes Characteristic Boundary Conditions (NSCBC) methodology and are based on the original formulation of [46] , incorporating the later improvements described in [47][48][49] that include source and transverse terms. Wall boundaries, where present, are treated as no-slip, isothermal, smooth solid surfaces and are implemented following the methodology described in [50] and [51] for non-porous, impermeable materials, such that the wall-normal mass flux of all chemical species is identically zero.

Spontaneous ignition of hydrogen reheat flames
In this section, we utilize one-and two-dimensional DNS calculations in order to investigate the ignition and stabilization characteristics of hydrogen flames at reheat conditions. The effect of the following parameters is studied: 1) the domain size ( L x ); 2) the fuel-oxidant equivalence ratio ( φ = f (Y k ) ); 3) the inlet velocity ( U in ); 4) the inlet temperature ( T u ); 5) the wall heat loss and flow confinement by walls (at y = 0 and y = L y ).

General features of the initial ignition process
The spontaneous ignition process is studied adopting an idealized (and simplified) representation of the reactive flow of interest. Figure 1 illustrates a sketch of the computational domain that is characterized by a (main) longitudinal dimension L x and, when present (e.g. in 2-D and 3-D calculations), by the transverse dimensions L y and L z (the latter is not shown and extends in the direction normal to the page, it is always a periodic/cyclic boundary). The reactant mixture consists of hydrogen and the vitiated oxidant stream originating from the first combustion stage. The latter is assumed to be a hot-gas generator that provides the products of hydrogen-air combustion at an equivalence ratio φ = 0 . 43 and reactants' temperature of T u = 773 K , mixed with additional air. Accordingly, the nominal target conditions of interest for the hydrogen reheat flame are defined, similar to [34] , as T u = 1100 K and φ = 0 . 35 (mass-fractions: 0.008 H 2 ; 0.183 O 2 ; 0.052 H 2 O; 0.757 N 2 ) resulting in nominal values for the ignition delay time, τ ig ∼ 0 . 15 ms , and adiabatic flame temperature, T ad ∼ 1800 K (from homogeneous reactor calculations). These are compatible with high efficiency and low emission in a typical gas-turbine combustion system (premixed).
In the present investigation, to find how the combustion may vary for conditions in the vicinity of the normal target condition, a parametric variation is introduced in the initial (fresh-reactant) temperature T u and composition Y k (i.e. equivalence ratio φ) of the reactant mixture that initially fills the entire domain at the start of each simulation. Entering into a continuous stream from the left (NSCBC inlet) boundary at x = 0 with a mean velocity U in , the reactant mixture is advected downstream until the initial spontaneous ignition occurs after a convective residence time t = t res ∼ τ ig , corresponding to spatial location x ∼ L ig,init ∼ U in · τ ig . The reacted, partially reacted or unreacted (depending on specific local conditions) gases exit the computational domain from the right (NSCBC outlet) boundary at x = L x . Note that, when ignition occurs, due to the way the simulations are initialized, the entire domain between  x = L ig,init and x = L x simultaneously ignites. This implies that the domain size affects the initial transient of the unsteady ignition process.
The initial transient phase of the spontaneous ignition process is described below for an exemplar (1-D DNS with the following parameters: φ = 0 . 12 , T u = 1100 K , L x = 10 cm and U in = 200 m/s ) as it is qualitatively similar for all cases investigated; however, quantitative differences emerge, depending on the specific conditions, and lead to different solutions that can remain unsteady or approach steady-state.
At a time t ∼ τ ig after the start of the simulation the spontaneous ignition process, in its initial phase that can be defined as "mildly explosive", leads to a sudden temperature increase throughout the downstream portion of the computational domain for L ig,init < x < L x . Locally, the temperature increase is accompanied by the simultaneous pressure increase and expansion of the gas mixture that is undergoing chemical reaction. This is graphically illustrated in Fig. 2 (a) by the dash-dot lines of the temperature (black) and pressure (red) profiles corresponding to the ignition time t = τ ig ∼ 0 . 2 ms . At later times during this initial explosive phase, the pressure wave generated by the spontaneous ignition process propagates downstream (exiting the computational domain through the NSCBC outlet) and upstream towards the fresh reactants and the domain inlet (exiting the computational domain through the NSCBC inlet). The latter, upstream-propagating pressure wave, however, as opposed to the downstream-propagating one, increases in amplitude as it moves towards the domain inlet ( Fig. 2 (a) covering the time interval t = 0 . 25 − 0 . 35 ms ) because it effectively represents an adverse pressure gradient for the approaching reactant flow (as an immaterial "piston" acting against it). The characteristic strength of the "piston effect" is directly proportional to the domain size (amount of reactants that ignite), equivalence ratio (temperature increase due to ignition) and reactant inlet flow velocity (steepness of adverse pressure gradient). The pressure increase that occurs in the fresh reactants also leads to a local temperature increase (up to a limiting value T u,max ) and to a shortening of the ignition time (down to a limiting value τ ig,min and a corresponding L ig,min ) that, in turn, causes the pro-gressive upstream displacement of the spontaneous ignition front (black lines in Fig. 2 (a)), strengthening the "piston" effect. The observed upstream combustion front displacement, although slower than the speed of sound (it lags the pressure wave), occurs against a mean flow U in ∼ 200 m/s . Accordingly, it is much larger than the "normal" flame deflagration velocity S L . Therefore, it can be concluded that the unsteady ignition phenomenon fits the criterion proposed by Zel'dovich for a spontaneous propagation regime.
Following the first, "explosive" phase of the unsteady spontaneous ignition, once the main pressure wave generated by the initial fluid expansion leaves the computational domain through the left boundary (NSCBC inlet) and the upstream propagation of the combustion front stops (due to τ ig (T u,max ) > t res ), the "piston" effect abruptly ends and this leads to a second, "relaxation" phase of the unsteady process. During the "relaxation" phase, the pressure decreases rapidly throughout the computational domain along with the temperature in the fresh reactants that decreases to its original value T u (set by the boundary condition). Spontaneous ignition is no longer sustained at all locations x < L ig,init and the combustion front is displaced downstream to its "natural" position where the solution reaches steady-state. This relaxation process is illustrated by the broken lines in Fig. 2 (b); please note that the solid lines represent the steady-state solution at the final (and much later) time t = 0 . 01 s .

Effects of the domain size and equivalence ratio
While the phenomenological picture described in the previous section qualitatively applies to all cases investigated, at least in respect to the first "explosive" phase, the transition to and the behavior of the second, "relaxation" phase is strongly dependent on the quantitative characteristics of the first phase, i.e. on the specific conditions of the simulated case. In the present section, the effects of domain size and fuel-oxidant equivalence ratio are discussed, while, in the next sections, an analysis of the effects of inlet-temperature variations is presented and the discussion extended to the case of confined flows with heat loss to the walls. Although the effect of the reactant inlet flow velocity on the igni- (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.) tion process is not explicitly illustrated and described here, it has been thoroughly investigated, and it is qualitatively similar to the effects of domain size and equivalence ratio described below, i.e. the strength of upstream propagating pressure waves increases for increasing inlet velocities. Figure 3 (a) illustrates, for eight different 1-D ignition cases, the time history of the maximum value of pressure recorded at each time step throughout the computational domain. A striking qualitative similarity is observed for all curves: a primary pressure peak due to the initial spontaneous ignition, followed by a secondary pressure peak due to the upstream-propagating pressure wave that exhibits a higher value compared with the former (the notable exceptions being the φ = 0 . 12 cases with the larger domains, L x = 20 cm and 30 cm ). The magnitude of the dual-peaked maximum-pressure time history shows a weak dependence on the domain size, i.e. larger domains give slightly higher pressure peak (compare the black, violet and pink lines in the figure), and a strong dependence on the equivalence ratio with the maximum pressure values increasing from P max ∼ 1 . 2 atm to ∼ 1 . 6 atm for an increase of the equivalence ratio from φ = 0 . 17 to 0.35 (blue, green orange and red lines). This is consistent with the increasing heat release with increasing equivalence ratio. It should be noted that for φ > 0 . 3 the displacement of the spontaneous ignition front, trailing the upstream-propagating pressure wave, rapidly reaches the one-dimensional domain inlet boundary and the NSCBC implementation is unable to handle the ensuing interaction, resulting in the simulations crashing at t ∼ 0 . 0 0 06 s . On the other hand, for φ < 0 . 3 , following the initial unsteady transient (consisting of the explosive and relaxation phases), all flames reach steady-state, with the spontaneous-ignition front positioned at the expected stabilization location L ig = U in · τ ig . Finally, before concluding the present section, it is important to mention that it is possible to stabilize hydrogen reheat flames at higher equivalence ratios (i.e. the target value φ = 0 . 35 ) by initializing the calculation with a reactant mixture at lower equivalence ratio i.e. φ = 0 . 17 and, after the initial transient is completed and the relaxation has occurred, increasing the amount of fuel introduced at the inlet boundary. The effect of this procedure is illustrated in Fig. 3 (b) that shows a gradual, smooth increase in the flame temperature following the increase in the hydrogen mass fraction imposed at the inlet boundary. Note that the maximum-pressure time history for this case, shown by the dashed blue line in Fig. 3 (a), is virtually indistinguishable from the case where no fuel increase is implemented, shown by the solid blue line. Figure 4 illustrates the typical de-stabilizing effect on the spontaneous ignition process of hydrogen reheat flames observed for values of the reactant mixture initial (inlet) temperature in the temperature range near crossover 980 K < T u < 1080 K (see Fig. 5 in [35] ). Although only results obtained for T u = 10 0 0 K are shown in Fig. 4 , the same trend is observed in all spontaneous-ignition tests conducted below a value of T u ∼ 1080 K . The first explosive phase of the spontaneous-ignition process is consistent with the description in the previous section for T u = 1100 K , see Fig. 4 (a). The relaxation phase, however, differs considerably. Here, a selfexcited instability of the flame emerges displacing the spontaneous ignition front back and forth in a nearly periodic fashion. A vast parametric study conducted in the framework of the present work for 980 K < T u < 1080 K and at atmospheric pressure conditions indicates that the amplitudes of the pressure and temperature fluctuations remain nearly constant (for the time duration of the DNS) or decrease for a reactant mixture equivalence ratio lower than Fig. 4 (c). Conversely, for equivalence ratios φ > 0 . 2 , the amplitudes of the pressure and temperature fluctuations exhibit non-monotonic and non-linear growth rates, see Fig. 4 (d). Note that flames are always stable for T u = 1100 K , as evidenced in Fig. 4 (c) and (d).

Response of hydrogen reheat flames to inlet temperature variations
At this point it is important to clarify the fact that the selfexcited flame instability observed for 980 K < T u < 1080 K is related not only to the unsteady initial ignition process. Even after stabilization of the spontaneous-ignition front is reached at the characteristic location L ig , the combustion process always can be destabilized by a reduction of the inlet temperature below T u ∼ 1080 K . This is illustrated in Fig. 5 (a) showing that, as soon as the inlet temperature (green curve) is reduced near the crossover temperature, a spatial oscillation of the spontaneous ignition front, with growing amplitude, occurs (in this case, causing a crash of the simulation). On the other hand, during an increase of the inlet temperature from T u = 1100 K to T u = 1150 K (away from crossover), no self-excited flame oscillation is observed, as shown in Fig. 5 (b). Additionally, the figure also illustrates that a later inlet temperature reduction from T u = 1150 K back to T u = 1100 K does not have a destabilizing effect, confirming that the observed self-excited combustion instability is not simply related to a generic reduction (in  itself) of the inlet temperature, but rather is induced by its decrease below 1080 K , towards the crossover value.

Temporal characterization of self-excited flame instabilities
The time scale of the observed periodic spatial oscillations of the flame, resulting from a self-excited intrinsic instability of the flame in the computational domain, is related to the convectiveacoustic feedback mechanism described by Williams, see p. 207 in [52] . Although the original reference is concerned with supersonic flows of reacting mixtures that involve shock waves with exothermic, finite-rate, temperature-sensitive chemistry occurring in the subsonic flow behind the shock, one particular and very prevalent mechanism described therein seems relevant to the present situation. A disturbance in the temperature of the incoming flow alters the auto-ignition time of the mixture, thereby perturbing the position at which ignition occurs. The perturbation in the ignition-front location produces a pressure pulse, which travels upstream acoustically, modifying the inlet temperature (just behind the shock, in the original configuration), which, in turn, further affects the autoignition time of the fluid element convected downstream. This then results in a self-sustained convective-acoustic feedback mechanism, the period of which is twice the sum of the time for a fluid element to be transported from the inlet to the ignition point (convective time t c ) and the time for an acoustic wave to travel from that point back to the inlet (acoustic time t a ). Twice because the oscillation involves a compression as the front moves upstream followed by a front-generated rarefaction as the front moves downstream.
For the specific example of the flame instability observed when T u = 10 0 0 K and φ = 0 . 12 and represented in Fig. 4 (a-c) of the convective-acoustic mechanism predicts that the oscillation period T W can be approximated as: resulting in a predicted T W ∼ 1 . 8 ms . The prediction agrees qualitatively with the observed oscillation period T f,obs ∼ 2 ms . The shortening of the oscillation period (i.e. frequency increase) that is visible in Fig. 4 (d) for the flame subjected to an increasing equivalence ratio is due to the movement of the mean location of the unstable flame during its oscillation cycle to spatial positions increasingly closer to the inlet boundary, such that the length scales that appear in the mechanism decrease with time. The heat release increase with time associated with the increasing equivalence ratio contributes to strengthening of the pressure waves which then are responsible for the observed increase of the amplitude over time. This is observed in the figure to be disrupted when the flame (temporarily) reaches the inlet twice.
In concluding the present section, it is important to highlight the fact that the self-excited instability of the flame presented above is also observed to take place in multi-dimensional configurations and in the presence of turbulence modulation. This is discussed in Section 3.4 for 2-D non-adiabatic configurations with quasi-realistic turbulent velocity fluctuations and in Section 4.3.2 for 3-D configurations with a realistic turbulent velocity field.

Effect of wall-confinement and heat-loss to the wall (non-adiabaticity)
The aim of the present section is to investigate whether the spontaneous ignition and flame stabilization processes described above for adiabatic conditions are significantly affected by the presence of non-adiabatic conditions. It is found that, in general, the observations presented previously are qualitatively valid also for the case of confined flows with wall heat loss, i.e. stable spontaneous ignition is achieved for φ ∼ 0 . 17 (initial equivalence ra-tio) and T u ∼ 1100 K while unstable ignition behaviour is observed above φ ∼ 0 . 18 and for T u < 1080 K .
Following the same wall boundary conditions implementation previously used in DNS studies performed with S3D [51,[53][54][55][56] , two isothermal, no-slip, smooth walls are placed opposite to each other separated by a distance of 1 . 5 cm in the transverse ydirection and kept at a fixed temperature T w = 750 K (i.e. lower than the fluid temperature) to form a 15 cm long straight channel flow where the spontaneous ignition process occurs. The relatively "cold" temperature of the isothermal channel walls (acting to confine the bulk flow of "hot" reactants) and the ensuing thermal boundary layers are intended to provide a simplified model representing the effect of wall-flushing by compressor air. Moreover, in order to represent the effect of turbulent convective heat transfer between the isothermal walls and the bulk flow within the intrinsic limitations of a two-dimensional configuration, a random flow field with a prescribed Passot-Pouquet energy spectrum (characterized by a rms velocity fluctuation of u = 25 m/s and an integral length scale of L T = 0 . 5 cm ) is superimposed onto the mean flow according to a well-established procedure [57] . The mean flow is described by a characteristic turbulent channel mean velocity profile with a centerline (inlet) velocity U c = U in = 200 m/s . The velocity fluctuations entering the domain from the inlet boundary, in conjunction with the wall heat loss, induce a temperature variance in the reactant mixture, acting to dissipate heat and radicals that are formed in the process that ultimately leads to a delay in spontaneous ignition. This is illustrated in Fig. 6 (a) by a comparison of the longitudinal temperature profile for the adiabatic onedimensional laminar case (black line) with the two-dimensional instantaneous (blue lines) and pointwise, time-averaged temperature profiles in the bulk flow (red lines). The wall heat loss and the temperature variance introduced by the relatively cold walls in the two-dimensional channel-flow configuration affects the nonadiabatic spontaneous ignition process that is displaced downstream by approximately 3 mm compared to the adiabatic process in the one-dimensional configuration. Note that, in the plots of Fig. 6 (a), the two 2 mm thick regions of the flow that are immediately adjacent to the isothermal channel walls are not considered. This finding confirms earlier observations about the role that temperature and compositional inhomogeneities play in delaying spontaneous ignition [58] .
Importantly, the occurrence of the characteristic self-excited intrinsic instability described in the previous section is also observed at non-adiabatic conditions for opportunely chosen values of T u and φ. A typical example is illustrated in Fig. 6 (b) showing the time history of the maximum pressure and temperature recorded within a L x = 15 cm long, wall-bounded duct crossed by a U c = 200 m/s mean flow for T u = 10 0 0 K and φ ∼ 0 . 185 . These conditions result in a nominal induction time τ ig ∼ 0 . 75 ms that is approximately equal to the channel residence time, t res = L x /U c . Accordingly, the nominal flame stabilization location is close to the downstream end of the duct. Note also that the first ignition event, taking place at t ig ∼ 1 . 3 ms , is delayed with respect to the nominal induction time for the mixture ( τ ig ∼ 0 . 75 ms ) due to the gradual increase of the equivalence ratio from φ = 0 at t = 0 to φ = 0 . 185 at t = 0 . 5 ms . The maximum temperature and pressure traces in Fig. 6 (b) clearly indicate the occurrence of a cyclic process characterized by distinctive time scales, while the sequence of plots presented in Figs. 7 and 8 illustrates the underlying spatial oscillation of the reaction front and the intermittent spontaneous ignition and extinction events. The latter are caused by the flame intermittently exiting the downstream end of the computational domain.
The two-dimensional simulations confirm the important role of compression heating, caused by the upstream-propagating pressure wave acting on the approaching reactant flow, during the first explosive phase of the spontaneous-ignition process. This is  front. This is quantitatively illustrated by the ratio of τ ig,loc /t res,loc ∼ 1 in the plots on the right-hand side of Fig. 7 , where the local value of the ignition time τ ig,loc is calculated by the analytic expression provided in [35] for hydrogen-air mixtures, while t res,loc is estimated as x/U c . The upstream displacement of the spontaneous ignition front halts once the ratio τ ig,loc /t res,loc (and its spatial gradient) becomes too large to be overcome by the effect of compression heating, see Fig. 7 (h). Figure 8 focuses on the subsequent relaxation phase, in which the reactant mixture temperature decreases (light green region immediately upstream of the spontaneous ignition front) because of the simultaneous local decrease in pressure, which in turn, is due to the abrupt interruption of the upstream displacement of the spontaneous ignition front. Auto-ignition cannot be supported any longer at such an advanced location due to the low(er) temperature of the approaching reactants. The reaction front is displaced downstream and, ultimately, flushed out of the channel (in the present configuration) before a new ignition cycle commences. This process results in intermittent spontaneous ignition cycles, alternating between ignition and upstream advancement followed by recession and extinction.
The characteristic time scales of the cyclic process observed in these two-dimensional DNS with wall heat loss are consistent with the convective-acoustic feedback mechanism of the combustion instability described in Section 3.3.2 for the one-dimensional adiabatic configurations. Approximating x f, 1 ∼ 5 cm , x f, 2 ∼ 15 cm , U in ∼ 200 m/s and c ∼ 650 m/s in Eq. (1) , a value of T W ∼ 1 . 44 ms is obtained which provides satisfactory agreement with the value T f,obs ∼ 1 . 6 ms observed in Fig. 6 (b) as the time period between each ignition and extinction event (the width of the red "bumps").

Background and rationale
Understanding the mean velocities of the reactants that are sufficient to stabilize a turbulent reaction front resulting from a hy-brid combustion mode of propagation and auto-ignition is important to the design of reheat combustion systems [59] . This is because mixed combustion modes transitioning from predominantly spontaneous ignition to predominantly flame propagation, are believed to limit engine operation with hydrogen-based fuels [12] . Therefore, using three-dimensional turbulent DNS, the main rationale for the present section is to obtain quantitative estimates of the ability of the turbulent reaction front of hydrogen reheat flames to balance the mean velocity of an approaching flow of reactants. Moreover, it is also shown that the self-excited intrinsic instabilities of the reaction front, described in Section 3 for oneand two-dimensional configurations, can occur in the modulating presence of a three-dimensional turbulent velocity field if specific conditions are met.
Although the present investigation builds upon the work by Savard et al. [33] , it substantially differs from their earlier work in several respects. Firstly, Savard imposes upon the mean flow and throughout the computational domain homogeneous isotropic turbulence supported by artificial (numerical) forcing of the largescale turbulent motions across the flame brush. This methodology was originally developed for isotropic, incompressible turbulence [60] and its applicability to combustion DNS remains controversial because of the implications of artificial forcing on the turbulencechemistry interaction dynamics, particularly on the burnt side of the flame which is not well understood as of yet. Secondly, Savard characterized the transition between spontaneous ignition and flame propagation utilizing relatively long computational domains that permit residence times of the order required for spontaneous ignition but, through the use of a low-Mach approximation, neglect compressibility effects and compression heating on the spontaneous ignition process. On the basis of the results described in Section 3 of the present study, there is compelling evidence suggesting that compressibility plays a significant role in hydrogen-reheat flame stabilization. The present DNS investigation uses relatively short computational domains and, placing the reaction front at a distance from the inlet boundary that theoretically can not support spontaneous ignition (initially), focuses on propagating turbulent reaction fronts and on their ability to balance the approaching upstream flow of the reactant mixture. Here, compressiblity effects are fully captured and able to detect eventual transition to mixed propagation/auto-ignition combustion modes due to local effects of compression heating, shortening of τ ig and subsequent early spontaneous ignition.

3-D DNS configuration
The approach chosen in the present DNS study employs an idealized reheat combustion configuration corresponding to a statistically planar turbulent premixed flame placed in unbounded turbulent flows of the target reactant mixture injected at the domain inlet boundary and characterized by auto-ignitive and adiabatic conditions (no heat loss). The three-dimensional computational domain has physical dimensions L x = 2 cm , L y = 1 cm and L z = 0 . 5 cm in the streamwise (non-homogeneous) x-direction, transverse (periodic/cyclic) y -and spanwise z -direction, respectively. The domain is discretized on a 800 × 400 × 200 Cartesian (uniform) mesh providing 25μm spatial resolution. The simulated three-dimensional turbulent hydrogen reheat flames represent the target conditions previously described in Section 3.1 ( T u = 1100 K , φ = 0 . 35 and T ad ∼ 1800 K at normal atmospheric pressure). For the present auto-ignitive conditions a unique reference laminar flame speed, S R , is defined as the inlet flow velocity at which the rate of change of the position of the flame from the inlet with inlet velocity is at a maximum [23] . This reference speed corresponds to the inlet flow velocity above which the steady-state laminar flame detaches from the domain inlet in a one-dimensional DNS configuration [23] . The reference speed and the associated thickness of an unstrained adiabatic laminar premixed flame in the limit of an unreacted upstream composition are S R ∼ 24 m/s and l R ∼ 1 . 35 mm (estimated by the maximum-thermal-gradient method). These values result in a reference chemical (flame) time scale, τ R = l R /S R = 5 . 63 e −2 ms . Notably, S R differs from the conventional "laminar-flame-speed" definition, S L , not applicable at auto-ignitive conditions, by accounting for the role of the residence time (owing to the relevance of ignition) on the flame speed.
In the main parametric investigation, the reactant flow is subjected to different levels of turbulence intensity ( Section 4.3 ) while in a second parametric sweep a different induction-time history of the flammable mixture ( Section 4.4 ) is investigated. The initially planar flames are subjected to different inlet turbulence intensities u = 3 (Case A), 10 (Case B), 15 (Case C) and 25 m/s (Case D). These turbulence intensities are specified as random 3-D flow fields with a prescribed Passot-Pouquet energy spectra, following a well-established procedure described in [57] . The velocity fluctuations are superimposed onto the mean flow that is advected into the domain from the upstream boundary with a velocity U in . The chosen conditions correspond to turbulent Reynolds numbers Re t = 22 , 75, 112 and 188, respectively. Assuming a size limit for the largest eddies in the flow equal to L T = 0 . 3 cm (well within the smallest transverse domain dimension, L z = 0 . 5 ), the corresponding longitudinal integral length scales, L 11 , lie between 0.11 and 0 . 31 cm . Table 1 summarizes the physical scales of fluid motion (al-ways referred to unburnt conditions) and of chemical reaction. It also provides estimates of the Karlovitz ( Ka ) and Damkøhler ( Da ) numbers for the simulated flames, placing them in the combustion regime diagram.
The turbulent flames are initialized in the 3-D computational domain through a progress-variable mapping from the corresponding laminar 1-D solution and placed at a distance x f = 0 . 5 cm from the domain inlet boundary. The progress variable C is a scalar parametrization of the reactive flow field, based for the present implementation on the hydrogen fuel mass fraction, that is equal to zero in the fresh reactants ( C = 0 ) and unity in the burnt products ( C = 1 ). An initial mean velocity U 0 ,in = 25 m/s is imposed at the domain inlet boundary and throughout the computational domain for all cases described in the following sections. For the duration of the simulation, the mean flow velocity imposed at the inlet boundary U in is adjusted at each time step such that the total amount of fuel that instantaneously enters the domain matches the volumetric fuel consumption rate of the deficient reactant, hydrogen. This procedure assures that the mean flow velocity U in ∼ U m is approximately equal to the displacement velocity S t of the turbulent flame reaction front, thereby ensuring that the latter remains within the computational domain at all times. This simple method is able to retain the mean flame position (approximately) in the vicinity of the initialization location x f , with only marginal upstream displacements. This is important for the following reasons: 1) the turbulent velocity fluctuations imposed at the inlet boundary are able to interact with the flame front before they are dissipated; and 2) the reactant-mixture residence time t res ∼ x f /U m between the inlet boundary and the flame position remains smaller than the ignition delay time, τ ig ∼ 0 . 15 ms , of the target mixture, thereby preventing a purely auto-ignition combustion regime from occurring.
The typical time evolution of the DNS solution is qualitatively illustrated in Fig. 9 , which provides a graphical representation of the turbulent (statistically planar) flame, at auto-ignitive conditions, as it responds to the approaching turbulent flow. The flow conditions represented in Fig. 9 , as an example, correspond to the turbulence intensity u = 15 m/s (Case C) and exhibit considerable wrinkling of the flame front (represented by the isosurface of temperature at T = 1500 K ) while, at the lowest turbulence intensity level u = 3 m/s (Case A), only a very mild wrinkling is observed (not shown). At the onset of the simulation, the initially flat surface marking the reaction-front location is wrinkled by the underlying turbulent flow field, and it rapidly accelerates, advancing in the upstream direction towards the inlet boundary. However, the procedure described above automatically adjusts the mean inlet velocity accordingly and ensures that the reaction front remains at a (statistically constant) distance from the inlet boundary plane.

Effect of the turbulence intensity
The main parametric investigation conducted quantifies the effect of the turbulence intensity of the approaching reactant flow on the turbulent reaction-front velocity at auto-ignitive conditions. Hydrogen flames at reheat conditions are characterized by values of the reference flame speed S R that are considerably larger than those found at "conventional" premixed-combustion conditions and, for the present configuration, S R ∼ 24 m/s .

Global analysis of the reaction-front velocity
The global reaction-front velocities, S t , estimated from the DNS cases, are presented in Fig. 10 in terms of non-scaled values. Scaled values are summarized in Fig. 13 (b) below. Time histories of the instantaneous (fluctuating) values of S t (red lines) for turbulent hydrogen reheat flames, subjected to inlet turbulence intensities u  Fig. 10 (a). The solution for a one-dimensional laminar configuration is also included for reference (black line). After an initial transient which occurs until time t ∼ 0 . 5 ms the reaction-front velocity, S t , relaxes towards its mean values of 32, 34, 36 and 42 m/s , respectively. The DNS are discontinued once a moving time window corresponding to 1 ms does not change more than 1% from the mean value of S t . Because the instantaeous fuel consumption rate in a turbulent flame is, inherently, a fluctuating quantity, in order to obtain a meaningful estimate for S t using this approach, the averaging must be conducted over a sufficiently long period. As discussed in [33] , it is unclear how long this period should be and, while previous DNS results suggest this period may be on the order of 10-100 τ T, 11 , for shorter averaging periods an uncertainty of approximately 5-20% on the calculated values of S t is proposed by Savard et al. (see Appendix B of [33] ).

Self-excited instability of turbulent hydrogen reheat flames
Although the rationale for performing the three-dimensional DNS principally concerns the estimation of the turbulent reactionfront velocity, in this section exploratory DNS are performed outside of the target (stable) conditions to investigate the occurrence of self-excited flame instabilities in the presence of a realistic, three-dimensional representation of the turbulent velocity field. As discussed in Section 3 , one-and two-dimensional DNS indicate that, at atmospheric pressure, self-excited flame instabilities typically arise at a reactants temperature T u ∼ 10 0 0 K due to the proximity to the crossover temperature. Therefore, an additional 3-D DNS, otherwise identical to Case B (i.e. u = 10 m/s ), is performed at T u = 10 0 0 K and atmospheric pressure (named Case B2). Moreover, three additional "variants" of Case B are performed for pressurized conditions, corresponding to P = 5 bar , at three reactants' temperature equal to T u = 10 0 0 K , 110 0 K and 1135 K (named B3, B4 and B5 respectively). The 3-D DNS configuration for the elevated pressure cases is formally identical to the atmospheric pressure reference Case B (i.e. same domain size, turbulence intensity etc) except for the finer mesh resolution ( 8 . 3 μm) that is required to resolve the smaller length scales of motion, diffusion and reaction of the reactive flow at elevated pressure. This implies, of course, a considerable increase in the computational cost of the pressurized cases and, therefore, these calculations are integrated for a shorter time interval ( 2 ms instead of 4 ms ).
The global turbulent reaction-front velocities, S t , estimated from these additional DNS datasets, "variants" of Case B, are presented in Fig. 10 (b) and compared to Case B, suggesting the following: 1. Onset of instability -At atmospheric pressure, the turbulent reaction-front velocity S t develops a distinct oscillation pattern with only minor spatial displacement (not shown) at T u = 10 0 0 K , Case B2 (black dashed line, delta symbols). This is shown by the almost sinusoidal pattern of the S t fluctuation and by a visibly larger oscillation amplitude compared to the reference Case B for which only small stochastic fluctuations can be observed (black solid line, gradient symbols). 2. Self-excited instability -At pressurized conditions, the turbulent reaction front develops a strong self-excited instability with significant upstream and downstream spatial displacement of the reaction front (see Supplementary material) at T u = 1135 K , Case B5 (red dashed line, diamond symbols). It is clear that the S t fluctuation exhibits higher order frequency harmonics. This is a sign of non-linear phenomena caused by the non-linear saturation of the flame at high amplitudes. 3. No instability -At pressurized conditions but lower reactants temperatures, T u = 10 0 0 K (Case B3) and 1100 K (Case B4), no self-excited flame instability is observed (solid and dash-dot lines, square and circle symbols, respectively). These two cases which are below the crossover temperature have sufficiently long induction times compared with the reactants residence times, and hence are predominantly controlled by deflagration rather than ignition-front propagation as in Case B5.
These observations are consistent with the findings of Section 3 and indicate that turbulent hydrogen reheat flames tend to exhibit self-excited instabilities for T u close to the temperature range near crossover ( T u < 1100 K at atmospheric pressure, T u > 1100 K at P = 5 bar , see Fig. 5 in [35] ). Interestingly, the insta- bility that arises in Case B2 appears to be weaker than those characterizing the "homologous" 1-D and 2-D configurations described in Sections 3.3.1 and 3.4 . This could be a consequence of turbulence modulation of the relevant processes that lead to the selfexcited flame instability but additional 3-D DNS would be required to clarify this point which is beyond the scope of the present study.
A plausible alternative explanation for the weak instability observed in Case B2 is also provided here. In the 3-D DNS configurations, the reaction front is initialized, by design, and stabilized, by the numerical procedure adopted, relatively close to the domain inlet boundary targeting a predominantly propagating flame, in accordance with the stated objectives of Section 4 . On the other hand, the 1-D and 2-D DNS configurations reported in Section 3 extend longitudinally 30 cm and 15 cm respectively (well beyond the 2 cm affordable in 3-D) and target predominantly autoigniting flames that are, by nature of the instability mechanism, more prone to develop stronger self-excited instabilities. Therefore, predominantly auto-igniting reaction fronts may develop stronger self-excited instabilities and higher turbulent reaction-front velocities compared to predominantly propagating ones. This hypothesis is supported by earlier findings about the prominent role of reaction (compared to diffusion) in the displacement speed budget for spontaneous ignition fronts [33,39] and by the present DNS results for the pressurized flames, shown in Fig. 10 (b). The occurrence of a strong self-excited instability is evident for the flame at the highest reactants temperature, T u = 1135 K (Case B5) that is both closest to the crossover temperature and characterized by the shortest induction time, τ ig ∼ 0 . 16 ms . This is comparable to the estimated mean residence time for the reactants t res ∼ S t /L x, f ∼ 0 . 2 ms (where L x, f indicates the estimated mean flame location) suggesting a predominantly auto-igniting reaction front.
Also in this three-dimensional configuration, the characteristic time scales observed in the strong instability of Case B5 are in good accordance with the convective-acoustic feedback mechanism of combustion instability described in Section 3.3.2 for onedimensional adiabatic configurations. Approximating x f, 1 ∼ 0 . 2 cm , x f, 2 ∼ 0 . 7 cm , U in ∼ 25 m/s and c ∼ 700 m/s in Eq. (1) , a value of T W ∼ 0 . 37 ms is obtained which provides satisfactory agreement with the value T f,obs ∼ 0 . 33 ms estimated in Fig. 6 (b) and observed more clearly from the maximum pressure trace (see Supplementary material).

Local analysis of the reaction front
The reason behind the observed increase in the displacement velocity of the reaction front for increasing turbulence intensity can be deduced by examining Figs. 11 and 12 . Increasingly strong wrinkling of the initially planar reaction front by the approaching turbulent flow causes an increase in flame surface area and, consequently, in the overall burning rate through a well-known mechanism that characterizes all turbulent premixed flames. However, the increase in flame wrinkling is also responsible for strong diffusion effects, focusing significant amounts of mobile, chemically crucial species (i.e. molecular hydrogen), radicals (i.e. the H-atom) and perhaps, to a lesser extent, enthalpy into pockets of unburnt reactants that, at reheat combustion conditions, are on the verge of undergoing spontaneous ignition. Enhanced concentrations of Hatom and enthalpy are found in these reactants' pockets located between increasingly curved portions of the flame front protruding towards its unburnt side where, instead, focusing of molecular hydrogen and local enrichment is typically observed.
On the burnt side of the flame brush the flame-stabilization location is largely unaffected by the turbulence intensity on average, as shown in Fig. 11 (a), which provides the streamwise profiles of the averaged temperature field (all profiles are essentially identical for x > 0 . 5 cm ). Conversely, on the unburnt side, the thickness of the preheat layer is strongly broadened by the increasing turbulence. Here, the averaging process used to construct the temperature profiles consists first of point-wise time-averaging (for the time period t = 1 − 4 ms ) followed by spatial averaging in the homogeneous y -and z -directions.
Crucially, the scatter plot of Fig. 11 (b) reveals significant effects of the turbulence intensity on the distribution in progressvariable space of the (very light) H-atom mass fraction and of the net production rate of HO 2 . At low values of the progress variable ( C < 0 . 1 ), beyond the main part of the preheat layer immediately upstream of the reaction zone, increasingly strong turbulence intensities result in higher H-atom mass fractions and higher positive values of the HO 2 net production rate (indicating enhanced production). This observation suggests increased diffusion towards C < 0 . 1 of H-atoms that quickly react to form HO 2 (via the elementary reaction H + O 2 + M < = > HO 2 + M ) centered around C = 0 . 01 as well as enhancing branching rates, especially at smaller values of C. Corresponding scatter plots of the H-atom reaction rate (not shown, see Supplementary material) exhibit fluctuations in its net production rate through chain branching that increase with increasing turbulence intensity, to such an extent that, at the highest intensity, there are even samples with net H-atom consumption around C values where the fluctuation of the HO 2 production rate is largest. Notably, the mean position of the peak heatrelease rate in progress-variable space, shown in Fig. 11 (c), is not affected by the turbulence intensity and is located at C ∼ 0 . 65 for all cases investigated. Analogous to the case of the mean flame position in physical space, illustrated by coincident temperature profiles for x > 0 . 5 cm in Fig. 11 (a), the mean temperature distribution in progress-variable space, shown in Fig. 11 (d), seems unaffected to any significant extent by the turbulence intensity for C < 0 . 1 , suggesting that heat conduction plays a minor role compared to the fast diffusion of the H-atom. Note also that, in the scatter plots of Fig. 11 , the sample envelope of the configurations characterized by higher turbulent intensity contains all those subject to lower intensities.
The effects of the turbulence-chemistry interaction process, mainly related to enhanced chain branching and HO 2 production in the fresh mixture by H-atom diffusion at high turbulence intensity, were inferred from the mean trends presented in Fig. 11 . These trends are further substantiated by the plots of Fig. 12 showing the instantaneous spatial patterns of key quantities in physical space. On the left column, the instantaneous isosurface of the progress variable at C = 0 . 01 is colored by the ratio D/R in Fig. 12 (a) and by the net reaction rate of HO 2 in Fig. 12 (c) where D represents the diffusion of sensible enthalpy (sum of Fourier and species transport term) and R the heat release due to combustion. The plots show a clear spatial correlation between the curved portions of the surface whose center of curvature is in the reactants with enhanced enthalpy due to diffusion (red-colored patches in Fig. 12 (a) and increased production of HO 2 radical (purple-colored patches in Fig. 12 (c)). The occurrence of spatially distributed, relatively high values of the HO 2 mass-fraction and of its production rate peaking in reactant pockets surrounded by the wrinkled flame front is illustrated in Fig. 12 (b) and (d), respectively, suggesting that (early) spontaneous ignition is occurring therein. Furthermore, the 2-D plots also suggest a causal relationship between H-atom diffusion (black vectors) into reactant pockets and the enhanced production of HO 2 . Figure 12 (e) confirms, on the curved portions of the C = 0 . 01 isosurface whose center of curvature is in the reactants, the occurrence of relatively high concentrations of H-atom (shown as red-colored elongated ridges on a mostly blue surface). These H atoms diffuse there from the hotter regions of the reaction front where they are produced, see Fig. 12 (f). Note the adoption, in Fig. 12 (e), of a vantage point on the opposite side of the flame compared to Fig. 12 (a) and (c) for better visualization of the spatial pattern of the H-atom mass fraction.
Summarizing, fast diffusion of a light radical species (the Hatom) and, to a lesser extent, of heat into reactant pockets encircled by the highly wrinkled flame are responsible for increased reactivity of the turbulent reaction front as a whole, leading to localized and intermittent spontaneous ignition that causes, ultimately, the observed increase in S t at higher turbulence intensities.

Effect of the mixture induction-time history
The present section describes an additional parametric investigation based on 3-D DNS with the objective to study the effects of the reactant-mixture induction-time history on the reaction-front displacement velocity S t , estimated using the same procedure described earlier. As opposed to Cases A-D presented in the previous section, where the reactant-mixture composition introduced at the domain inlet contains only the major species (H 2 , O 2 , N 2 , H 2 O) and effectively represents an unreacted state, the DNS calculations presented in this section impose a partially reacted inlet-mixture composition and temperature at the domain boundary. Two compositional variations are considered for this partially reacted mixture, and these are extracted from a 0-D homogeneous-reactor calculation that is initialized with the target conditions defined in Section 3.1 and that result in the "standard" ignition time delay, τ ig ∼ 0 . 15 ms . The first mixture composition and temperature are extracted from the 0-D homogeneous reactor calculation after a time, t = 0 . 055 ms ∼ 3 / 8 τ ig , and the resulting partially reacted state is denoted here as "Advanced Ignition1" (AdvIgn1), while the second mixture composition and temperature are extracted after a time, t = 0 . 11 ms ∼ 3 / 4 τ ig , and the corresponding state is denoted as "Advanced Ignition2" (AdvIgn2). The rationale for this investigation lies in the importance of the "history" of different fluid parcels, or streamlines, in defining an effective flame front burning rate which actually results from the combination of varying local conditions. For an application within the industrial burner geometry, indeed characterized by local variations of the fuel mixture fraction and temperature, the mixing section cannot be considered as one single element at a single temperature and composition, but, based on the outcome of this study, the detailed modelling of the mixing section can be improved, considering the actual impact on the flame front of fluid parcels with different histories in the definition of an effective burning rate.
Results from the three-dimensional DNS are presented in Fig. 13 and indicate a considerable departure of the reaction-front displacement velocity S t when partially reacted mixtures are introduced at the domain inlet boundary (orange/red lines and symbols) compared to the reference case of an unreacted state for the inlet mixture (black lines and symbols). For the "advanced ignition" states of the mixture compositions and chemical reactions a clear increase in the overall fuel consumption rate within the reaction front occurs, leading to a greater displacement velocity of the front against the approaching turbulent flow. In proportion, the mixture at the more advanced ignition state (AvdIgn2) has a much larger impact on S t compared to the mixture at the less advanced ignition state (AdvIgn1) as shown in Fig. 13 (a). Figure 13 (b) compares the estimates of S t , obtained from the present calculations and normalized by S R (rectangles), versus a simple scaling law (solid line). The figure also illustrates that the present S t estimates are located well below the empirical law (dashed line) suggested by Savard et al. [33] to distinguish between combustion regimes dominated by ignition fronts (above the dashed line) or deflagrations (below it). Here, caution should be exercised in the interpretation of the S t estimates for the more advanced ignition state (AdvIgn2) state because: 1) the slight slope in the red-colored curves shown in Fig. 13 (a) suggest that convergence to a steady mean value for the S t estimate is not yet reached; 2) inspection of the instantaneous flame position in the 3-D DNS dataset (not shown) indicates that the turbulent flame brush is located very close to the domain inlet boundary and (undesirable) interaction of the reaction front with the boundary condition formulation cannot be ruled out for AdvIgn2. Finally, an interesting observation is also related to the effect of the turbulence intensity, that is still present, and that contributes to a further velocity increase, in addition to that caused by the "pre-cooking" of the reactant mixture introduced at the domain inlet boundary. This supports the interpretation of the turbulence effect on the flame front reactivity explained in the previous subsection.

Conclusions and further work
A series of Direct Numerical Simulations were performed to study several aspects of hydrogen-air combustion at auto-ignitive conditions (e.g. reheat combustion). The physical characteristics of the unsteady spontaneous ignition process in hydrogen reheat flames have been investigated in detail for adiabatic laminar conditions (in 1-D configurations) and for turbulent conditions, with and without wall heat-loss, (in 2-D/3-D configurations). The initial and boundary conditions that lead to a stable or unstable spontaneous ignition process have been identified. In addition to the reactant-mixture equivalence ratio φ, characterized by a threshold between 0.1 and 0.3 for a sufficiently weak (and controllable) ig-nition, the stability of hydrogen reheat flames is strongly affected by the premixture temperature. Near the crossover temperature, unsteady spontaneous ignition fronts always exist. An analytic expression is derived that satisfactorily approximates the time scales characterizing the unsteady motions of the reaction front observed in one-, two-and three-dimensional DNS calculations.
Furthermore, 3-D DNS of turbulent hydrogen-air statistically planar flames have been performed in order to obtain quantitative estimates of the velocity S t that the reaction front is able to sustain against the approaching turbulent flow at auto-ignitive conditions. The 3-D DNS have principally covered a range of turbulence intensities and mixture induction-time history for reactants at a temperature T u = 1100 K and at atmospheric pressure. Several additional exploratory DNS for different values of T u and at elevated pressure were conducted to confirm the occurrence of self-excited flame instabilities in turbulent hydrogen-air flames at auto-ignitive conditions. Moreover, results indicate that localized diffusion of fast diffusing H atoms and, to a lesser extent, heat from the wrinkled flame surface into reactant pockets plays a key role in increasing the reactivity of the entire reaction front as a function of increasing turbulence intensities. Values of S t between 30 and 60 m/s are observed at the conditions investigated implying that, in the gas turbine sequential combustion-system design, flow velocities well in excess of 60 m/s across the mixing section of the burner, spanning the bulk region of the flow, all the way to the boundary-layer regions, must be implemented in order to ensure robust hydrogen reheat flame stabilization and to avoid upstream displacement of the reaction front (i.e. flashback). The present findings could provide a plausible explanation for the recent observation, in full-size full-load experiments on the GT36 combustion system [12] , of a transition from autoignition stabilization to a hybrid flame-propagation mode that occurs for hydrogen concentrations in the fuel exceeding 70% (on a volume basis, the rest being natural gas). Proceeding beyond the atmospheric pressure conditions investigated here, future work will investigate the effect of higher pressure levels on the ignition and propagation characteristics of hydrogen reheat flames and of turbulent modulation on the self-excited reaction-front instability.

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