Edwin R C Baynes
email: edwin.baynes@univ-rennes1.fr

Dimitri Lague

Jean-Jacques Kermarrec

Supercritical river terraces generated by hydraulic and geomorphic interactions

The alternating cycle of valley widening through lateral erosion ('strath planation') and valley narrowing through vertical incision into bedrock ('strath terrace abandonment') due to variations in sediment supply (Q¬s) relative to river transport capacity (Qsc) is a common feature in many mountainous environments, yet our understanding of the mechanics of the processes that drive this landscape change remains poorly quantified. Here, we used an experimental and numerical study to identify the geomorphic and hydraulic controls driving the response of mixed bedrock-alluvial rivers to variable sediment supply, water discharge and tectonic tilting. The experimental channels exhibit a multi-stage response of channel narrowing, stripping of the alluvial cover in a downstream migrating incision wave followed by destabilization of the bed and development of a single vertical step in the bed profile ('knickpoint') when the hydraulic conditions are supercritical. In our experiments headward erosion by knickpoints is the most efficient process of strath terrace abandonment, contributing the majority of the total vertical incision in a short period of time. We show experimentally, that knickpoint developing under supercritical flow conditions, driving the rapid response of fluvial systems to upstream perturbations in Qs/Qsc despite no base-level fall. This has implications for the understanding of distributions of strath terrace ages, the inference of baselevel variations from knickpoint propagation, and how landscapes respond to climatic or tectonic perturbations.

INTRODUCTION

Past cycles of global climatic change have had consequences for the hydrological regime and the pattern of sediment production and transport within landscapes [START_REF] Schumm | Implications of Complex Response of Drainage Systems for Quaternary Alluvial Stratigraphy[END_REF][START_REF] Fuller | Field evidence for climatedriven changes in sediment supply leading to strath terrace formation[END_REF][START_REF] Wobus | Does climate change create distinctive patterns of landscape incision?[END_REF]. The balance between the flux of sediment entering a mixed bedrock-alluvial river (Qs) and the capacity of the river to transport that sediment (Qsc) is important in setting the landscape morphology, such as the generation and preservation of flights of strath terraces above the active channel (e.g., [START_REF] Pazzaglia | Bedrock fluvial incision and longitudinal profile development over geologic time scales determined by fluvial terraces[END_REF][START_REF] Wegmann | Holocene strath terraces, climate change, and active tectonics: The Clearwater River basin, Olympic Peninsula[END_REF].

Several mechanisms have been proposed for the generation of strath 'terraces', including climate driven oscillations in Qs/Qsc resulting in valley widening (high Qs/Qsc) or vertical incision (low Qs/Qsc) [START_REF] Bull | Stream-terrace genesis: Implications for soil development: Geomorphology[END_REF][START_REF] Hancock | Numerical modelling of fluvial strath-terrace formation in response to oscillating climate[END_REF]. The presence of strath terraces within a fluvial system have also been shown to result from headward incision by 'base-level knickpoints' initiated by fault throw or tectonically driven base level fall [START_REF] Finnegan | Interpretation and downstream correlation of bedrock river terrace treads created from propagating knickpoints[END_REF][START_REF] Finnegan | Sediment supply, base level, braiding, and bedrock river terrace formation: Arroyo Seco[END_REF], unsteady lateral river erosion [START_REF] Limaye | Numerical model predictions of autogenic fluvial terraces and comparison to climate change expectations[END_REF], bedrock structure [START_REF] Wohl | The effect of bedrock jointing on the formation of straths in the Cache la Poudre River drainage, Colorado Front Range[END_REF], the cutoff and migration of river meander bends [START_REF] Finnegan | Episodic bedrock strath terrace formation due to meander migration and cutoff[END_REF], sea level highstands [START_REF] Merritts | Long river profiles, tectonism, and eustasy: A guide to interpreting fluvial terraces[END_REF], river capture [START_REF] Stokes | Quantification of river-capture-induced base-level changes and landscape development[END_REF]. The influence of tectonics in generating strath terraces is usually limited to the consideration to the initiation of headward migrating incision waves (e.g., Finnegan 2013), while the impact of regional landscape tilting is rarely considered. A clear identification of the mechanism of strath terrace formation is important for the interpretation of the presence and pattern of strath terraces within a landscape to date past climate changes (e.g., [START_REF] Fuller | Field evidence for climatedriven changes in sediment supply leading to strath terrace formation[END_REF][START_REF] Stokes | Controls on dryland mountain landscape development along the NW Saharan desert margin: Insights from Quaternary river terrace sequences(Dades River, south-central High Atlas, Morocco)[END_REF], estimate fluvial incision rates (e.g., [START_REF] Lavé | Fluvial incision and tectonic uplift across the Himalayas of central Nepal[END_REF], or to help calibrate empirical models of fluvial erosion [START_REF] Finnegan | Interpretation and downstream correlation of bedrock river terrace treads created from propagating knickpoints[END_REF]. However, a universal model of strath terrace formation is currently limited by an incomplete understanding of the mechanics of vertical and lateral erosion in bedrock and difficulties in identifying the erosional pattern during the response of channel processes to perturbations in Qs/Qsc from preserved strath terrace age distributions.

Here, we present a physical modeling study that explores the mechanisms of channel response and resulting terrace geometry to three scenarios of perturbations to Qs/Qsc. The first two experiments explore the response of a channel following a sudden perturbation to Qs/Qsc from the same initial conditions while the third explores the impact of a gradual decrease in Qs/Qsc, through progressive tilting of the experiment. The experiments identify patterns of erosional response including the resulting preserved terrace morphology, aiding the interpretation of terrace age distributions in natural landscapes for purposes such as the reconstruction of climate or tectonic histories (e.g., [START_REF] Fuller | Field evidence for climatedriven changes in sediment supply leading to strath terrace formation[END_REF][START_REF] Stokes | Controls on dryland mountain landscape development along the NW Saharan desert margin: Insights from Quaternary river terrace sequences(Dades River, south-central High Atlas, Morocco)[END_REF]

EXPERIMENTAL SETUP

Experiments have been used extensively to explore the dynamics of bedrock channels, including abrasion processes (e.g., [START_REF] Johnson | Feedbacks between erosion and sediment transport in experimental bedrock channels: Earth Surface Processes and Landforms[END_REF], the dynamics of bedrock plucking (e.g., [START_REF] Dubinski | Relationships between block quaryying, bed shear stress, and stream power: A physical model of block quarrying of a jointed bedrock channel[END_REF], the mechanics of waterfall erosion in homogenous bedrock (e.g., Baynes et al., 2018), and landscape stability (e.g., [START_REF] Hasbargen | Landscape instability in an experimental drainage basin[END_REF]. We used the 80 cm long, 30 cm wide and 20 cm deep Bedrock River Experimental Incision Tank at the Université de Rennes (Fig. 1; S1), with an initial configuration of a thin (2 cm) layer of well-sorted sand (grain size = 250 µm) overlying homogenous silica paste (grain size = 45 µm). The cohesive silica paste material is commonly used to represent bedrock and is eroded by hydraulic shear under clear water flow in the absence of sediment (e.g., [START_REF] Turowski | Dynamics and steady-state geometry of an experimental channel incising cohesive material[END_REF][START_REF] Turowski | Experimental channel response to tectonic uplift[END_REF][START_REF] Turowski | Controls on bedrock channel morphology: experimental and theoretical investigations and comparison with natural channels[END_REF] and the experimental channels in this study exhibit a similar hydraulic scaling of width and slope to natural bedrock rivers (Baynes et al., 2018). No direct scaling between our experiments and a particular target natural system is sought but the behavior of the experiment is similar to that observed in natural systems, allowing the dynamics of processes and their controls to be explored quantitatively (Hooke, 1968; SI Section 1).

Constant input water discharge, set within 2% by a regulated pump on the inflow pipe, was fed into a reservoir upstream of the inlet, and sediment was added using a Gericke infinite screw feeder system. The channel outlet elevation was fixed, although it was free to migrate across the width of the box, and sediment was collected in a reservoir beyond the outlet to avoid deposition in the channel as a result of the outlet. The same initial conditions for the first two experiments were Q = 1.5 l/min, Qs = 13.33 g/min, slope = 4.5°, selected after previous experiments for the best combination of parameters to generate an alluvial river at equilibrium conditions inset in the sand material. Channel slope, width and depth were free to evolve during the course of the experiment through erosion (vertical and lateral) and deposition. After 37 minutes, once the channel had reached a constant width, one parameter was altered while all others kept constant to represent a sudden shift in Qs/Qsc; experiment 1: Qs reduced to 0 g/min, experiment 2: Q doubled to 3 l/min. In the third experiment, Qs and Q were maintained constant throughout but the slope was increased from an initial slope of 2.5° at a continuous rate of 1° per hour for 4.7 hours, gradually decreasing Qs/Qsc. A Leica ScanStation 2 Terrestrial Laser Scanner (green laser footprint = 3 mm) was programmed to collect point clouds of the topography every 2 min. The topography was then gridded with a pixel size of 2 mm, and combined with the hydrodynamic model Floodos (Davy et al., 2017), to evaluate the hydraulic conditions within the channel (e.g., Froude number) and to quantify the evolution of the channel geometry (e.g. width and depth) during the experiments (SI Section 2).

EXPERIMENTAL RESULTS AND DISCUSSION

Changing Qs/Qsc had a significant, rapid and complex impact on the initial stable state (channel width: 150 mm) in the first two experiments. Following the removal of input Qs in experiment 1, an initial narrowing of the channel coincided with a downstream migrating incision wave initiated at the inlet that removed the alluvial cover (Fig. 2A), exposing the bedrock surface to erosion. The removal of input Qs led to Qs/Qsc = 0, providing capacity for the flow to entrain and strip the sediment stored in the alluvial bed. Further narrowing of the channel occurred until the bed destabilized into a series of steps (spaced at 25-30 mm with height of 3-5 mm, at t = 83 min) in the bedrock surface that evolved into a single vertical step (15-20 mm in height at t = 111 min) that migrated upstream as a knickpoint, depositing material downstream of the plunge pool (Fig. 2A). The same multi-stage process of downstream stripping of alluvial cover followed by knickpoint development and headward erosion also occurred in experiment 2 where Q was doubled (SI Section 3). The development of steps before the generation of a single knickpoint were not observed in the topographic data from experiment 2 as the channel evolved from downstream stripping to headward erosion within the time period between successive laser scans (t = 120 s; Fig. S5B). Despite this, the knickpoint formed mid-channel in the absence of any base level fall or tectonic uplift.

The initial stages of experiment 3 (until t ~100 min) were characterized by an alluvial fan-like regime near the inlet with cone-shaped deposits frequently reworked by high lateral channel mobility (Fig. 2B) while the slope increased (Fig. S5C). As Qs/Qsc decreased, the flow focused into a narrower single thread, triggering the onset of a downstream incision wave of the alluvial cover. Subsequently, a knickpoint developed in the bedrock (t = ~208 min) which retreated upstream in the same manner as the other experiments (Fig. 2B; SI Section 3). In each experiment, the stripping of alluvial cover by the downstream incision wave coincided with a significant narrowing of the channel inset in the alluvial cover, with the narrowest width occurring in bedrock when a knickpoint had formed (Fig. 2B). The channel narrowing associated with the incision of both the alluvial cover and the bedrock led to the abandonment of the formerly active bed at higher elevations, representing the generation of strath terraces seen in natural environments. However, the rates of vertical incision and terrace abandonment were not spatially or temporally constant during the experiments, with the highest local rates of vertical incision associated with the downstream stripping of sediment and the headward migration of the knickpoints through bedrock (Fig. 2B; SI Section 4). A similar contribution to total vertical incision by knickpoint migration was identified in Hawai'i [START_REF] Mackey | Knickpoint formation, rapid propagation, and landscape response following coastal cliff retreat at the last interglacial sea-level highstand: Kaua'i, Hawai'i[END_REF].

The Floodos hydrodynamic model output (SI Section 5; Fig. S10) shows that the sequence of steps, and subsequently the knickpoint, developed in experiment 1 where flow had a Froude number > 1. The Froude number, a dimensionless number defined as the ratio of the inertia forces to gravitational forces, characterizes the transition from subcritical to supercritical flow when Fr > 1. During supercritical flow, interactions between the flow and the bed are inherently unstable [START_REF] Yokokwa | Cyclic steps by bedrock incision[END_REF]Izumi et al., 2017;[START_REF] Scheingross | A mechanistic model of waterfall plunge-pool erosion into bedrock[END_REF], and the development of steps in the bedrock surface occur during our experiments when the flow conditions exceed the threshold Froude number for supercritical flow. We therefore refer to these features as 'supercritical knickpoints' to distinguish them from 'base-level knickpoints'. Subsequently, we refer to the terraces generated by the passage of supercritical knickpoints as 'supercritical terraces'.

IMPLICATIONS FOR THE UNDERSTANDING OF NATURAL LANDSCAPES

Our experiments provide direct evidence of the hydraulic processes that lead to strath terrace formation following perturbations to Qs/Qsc. [START_REF] Wobus | Does climate change create distinctive patterns of landscape incision?[END_REF] show that climate change induced incision was associated with a downstream incision wave due to changes in the relative sediment supply to transport capacity, while base-level induced (e.g, river capture, tectonics) incision was characterized by an upstream migrating incision wave.

Crucially, our experiments show that channels in homogenous substrate respond to climatic perturbations through a complex multi-stage incision mechanism, containing a downstream (Fig. 3A) and upstream (Fig. 3C) wave component that both generate terraces, in the absence of any base level fall.

The multi-stage nature of the incision (Fig. 3) has important implications for studies that use the age of terraces to reconstruct climatic changes (e.g., [START_REF] Fuller | Field evidence for climatedriven changes in sediment supply leading to strath terrace formation[END_REF] or rates of fluvial incision and tectonic uplift (e.g., [START_REF] Lavé | Fluvial incision and tectonic uplift across the Himalayas of central Nepal[END_REF], due to the time transgressive nature of the terraces generated by this response (Zaprowski et al., 2001). During the downstream incision wave, the upstream parts of the alluvium are abandoned first, while the opposite is true during the upstream incision wave, leading to spatially variable terrace ages (Fig. 3D). If it is incorrectly assumed that the entire length of a terrace was abandoned at the same time, the average vertical incision rate calculated from a single terrace abandonment age would be misleading due to the failure to identify the erosion history of an incision wave and the spatially and temporally non-uniform nature of the river incision [START_REF] Finnegan | A signature of transience in bedrock river incision rates over timescales of 104-107 years[END_REF]. Strategic sampling of high-resolution point ages from terraces may differentiate between strath terrace formation by vertical incision or knickpoint retreat (e.g., samples collected longitudinally from the same terrace surface; [START_REF] Baynes | Erosion during extreme flood events dominates Holocene canyon evolution in northeast Iceland[END_REF], allowing accurate interpretation of the mechanisms and triggering perturbations that generated them.

Knickpoints in a channel are traditionally interpreted as markers of either local lithological boundaries (e.g., [START_REF] Haviv | Evolution of vertical knickpoints (waterfalls) with resistant caprock: Insights from numerical modelling[END_REF], a history of base level fall or tectonic uplift rate variability (e.g., [START_REF] Finnegan | Interpretation and downstream correlation of bedrock river terrace treads created from propagating knickpoints[END_REF], or threshold drainage areas for channel incision [START_REF] Crosby | Knickpoint intiation and distribution within fluvial networks: 236 waterfalls in the Waipaoa River, North Island[END_REF]. However, our experiments indicate that the presence of a knickpoint in a channel with homogenous bedrock and constant baselevel is not necessarily diagnostic of any of these sources as knickpoints are generated following perturbations to Qs/Qsc (Fig. 3).

Under normal conditions, flow in large natural alluvial rivers where strath terraces are commonly found is subcritical, but above-threshold supercritical flow conditions (Fr > 1) can occur due to a combination of driving factors including channel narrowing with variable width and increased discharge (SI Section 5). Base-level terraces are likely to have the same or a similar slope to the active channel bed through parallel knickpoint retreat [START_REF] Finnegan | Interpretation and downstream correlation of bedrock river terrace treads created from propagating knickpoints[END_REF], whereas terraces generated during the multi-stage response identified here are steeper than the active channel (Fig. 2, 3) due to new equilibrium conditions with lower sediment supply, higher discharge or continued tilting. An additional diagnostic link could be the highest terrace having a decreasing abandonment age with distance downstream (Fig. 3D) that does not occur during terrace formation driven by base-level fall. Thus, the morphology of terraces and spatial patterns of abandonment ages can be diagnostic of the processes and forcing perturbations that led to terrace generation and should be considered carefully when analyzing topography for the purpose of understanding and reconstructing past tectonic or climatic histories.

CONCLUSIONS

Landscapes are able to respond rapidly following climatic or anthropogenic perturbations that switch the erosion processes from one regime to another. The transition from lateral erosion and strath planation to channel narrowing, incision and strath terrace abandonment occurs through a multi-stage process of downstream then upstream incision waves. The headward passage of knickpoints can be initiated even in the absence of base level fall, due to the interaction between the channel bed and supercritical flow conditions. The passage of supercritical knickpoints provide an efficient process for the abandonment of strath terraces, and should be considered when using the presence of strath terraces within landscapes to infer past fluvial incision rates or climatic changes, which may be more spatially and temporally variable than previously thought. Experimental modelling studies of landscape and fluvial dynamics are typically performed using one of two approaches. The first directly scales the experimental prototype to the natural system by ensuring that non-dimensional hydraulic parameters such as the Froude number are scaled strictly and the physical processes (e.g., sediment transport) are the same in both the experimental and natural systems. The second approach relaxes this strict scaling, and seeks a 'similarity of process' between the experiment and the natural channel (Hooke, 1968) whereby the behaviour of the processes within the experiment system are qualitatively similar and 'reasonably effective' at replicating the behaviour of natural systems despite large differences in the governing dimensionless numbers [START_REF] Paola | The ``unreasonable Effectiveness'' of Stratigraphic and Geomorphic Experiments[END_REF]. This second approach has the benefit of allowing temporal and spatial timescales to be relaxed, as materials such as silica cement that is eroded by clear water flow can be used to replicate bedrock, and experiments exploring bedrock erosion processes can be performed relatively quickly. It should be noted that the results from these experiments cannot be directly scaled up to any particular natural system. Additionally, the experimental channel is always 'active', and the experiments do not include the periods of time between flood events when the channel is not being eroded, creating a difficulty when exploring relative timescales of geomorphic response between the experiments and natural rivers. However, processes such as flow acceleration above the knickpoint lip, cyclic step development, plunge pool erosion, undercutting of the knickpoint face and channel banks, cantilever failure, erosion and transport of cohesive material by hydraulic shear are present in the experiments. The typical width-depth ratio of the experimental channels at the start of the experiments (alluvial channel conditions) is 50-80, similar to that of natural gravel-bedded rivers [START_REF] Finnegan | Controls on the channel width of rivers: Implications for modelling fluvial incision of bedrock[END_REF], and the scaling of the channel width with discharge observed in natural rivers is also reproduced in the experiment. Therefore, we can be confident that the experiments are 'reasonably effective' in replicating natural processes, and the response of the channels to the perturbations (sediment supply, discharge or slope) in the experiments can help inform the response of natural channels to similar forcing.

FIGURES

Figure S1 shows the setup used during our experiments in the Bedrock River Experimental Incision Tank at the Université de Rennes. The initial alluvial cover was composed of wellsorted sand (grain size = 250 µm), and the homogenous cohesive silica paste used as a bedrock simulant was composed of 61.5% angular silica, 20.5% spherical beads and 18% water. The silica paste was mixed in a cement mixer before being transferred into the box flume and homogenised using a high frequency vibrating needle to re-liquefy the paste (see Baynes et al., 2018 for further details related to the properties of the silica paste and the preparation procedure). Flow in the channel is typically laminar with low Reynolds numbers (< 1000), and is input into the channel through a reservoir and the discharge magnitude controlled using a pump.

Fig. S1: Photograph of the Bedrock River Experimental Incision Tank at Université de

Rennes. The blue arrow shows the flow direction within the channel.

SI Section 2: Floodos parameters

Floodos is a numerical precipiton-based hydrodynamic model (see Davy et al., 2017 for full technical description of the model), and was used here to assess the hydraulic conditions during the experiments. The input topography for Floodos was produced directly by gridding the point clouds, collected from the terrestrial laser scanner at regular time intervals during the experiments, with a pixel size of 2 mm. Floodos was run using the laminar flow option consistent with the relatively low Reynolds number of the experiment (< 1000) across most of the channel under alluvial conditions (Fig. S2).

Fig. S2: (A) Digital Elevation Model from Experiment 1, when the channel is under fully alluvial conditions (t = 33 minutes). (B) Modelled water depth mask using Floodos

hydrodynamic model (Davy et al., 2017). (C) Calculated Reynolds number for each pixel of the flow, using the Floodos model output.

Under the laminar flow conditions, the local vertically averaged velocity U is approximated by: and should be approximately C~ρg/3μ where ρ is water density and g is gravitational acceleration. This predicts that at 10°c, C ~ 2.5 x 10 6 m -1 .s -1 . Fig. S2 shows the impact of different values of C on the predicted extent of water flowing over a flat surface compared to manual measurements of the wetted area obtained from a photograph taken from above the experiment (Fig. S3A). The simulations of flow extent are relatively similar for all modelled values of C, although lower values of C predict a greater level of overflow on the edges in the downstream part. We therefore used the value of C ( 2.5 x 10 6 m -1 .s -1 ) for the simulations presented here, due to good consistency with the manual measurements and the theoretical prediction of the friction coefficient based on the water viscosity for flow over a flat surface (Fig. S3) and for channelised flow (Fig. S4).

𝑈 = 𝐶𝐻 2 𝑆, (

SI Section 3: Evolution of channel and terrace geometry during experiments

This SI section provides the complete experimental data for the evolution of the channel geometry for each of the three experiments. The long profiles were extracted by finding the point of minimum elevation for each row of the grid (Fig. S5). Terraces were defined from the DEM as flat surfaces (elevation within 1mm) and at least 2 cm wide. The channel crosssections (Fig. S6) were extracted from the data at the location of the grid that was 25% of the distance from the inlet to the outlet of the channel. and Izumi, 2000). In this section, we used the Floodos hydrodynamic model output to calculate the Froude number for all locations within the channel (Fig. S8) for selected times of experiment 1 (t = 81 min, when bed instabilities were developing; Fig. S8A, and t = 106 min, when a single knickpoint had formed; Fig. S8B).

Cyclic steps form under supercritical flow conditions due to an interaction between the channel bed and the flow above it (Izumi et al., 2017), developing a series of steps with transcritical flow that alternates between supercritical and subcritical flow with hydraulic jumps (Fig. S8A). After a single knickpoint formed, a pronounced change in the Froude number and a hydraulic jump is present at the knickpoint lip, with subcritical flow conditions present in the plunge pool downstream (Fig. S8B). The absolute values of Fr are possibly slightly overestimated due to the use of the laminar flow option in Floodos, and some uncertainty related to whether the flow is purely laminar or transitional (Fig. S2). However, the pattern of Fr variation is consistent with cyclic step development as observed and predicted theoretically under supercritical flow conditions (Parker and Izumi, 2000;Izumi et al., 2017).

Supercritical flow in natural rivers

Here, we performed a simple test to consider whether supercritical flow conditions can be sustained in natural rivers under steady, uniform flow, by exploring the impact of varying channel width and discharge on the Froude number, assuming a rectangular channel with a constant slope.

The flow velocity can be estimated using the Manning's equation (U = n -1 .H 2/3 . S 0.5 ), where n is the Manning's roughness coefficient, H is the flow depth, and S is the channel slope. The flow velocity is also given by U = Q / WH, where Q is the discharge and W is the channel Eq (S3) shows that in steady uniform flow conditions, the Froude number is mostly sensitive to slope and roughness, with a weak dependency on width and discharge.

According to [START_REF] Parker | Physical basis for quasiuniversal relations describing bankfull hydraulic geometry of single-thread gravel bed rivers[END_REF], the channel width and channel slope of a gravel bed river at bankfull discharge conditions can be calculated using the following equations:

𝑊 𝑏𝑓 = Where Qbf is the bankfull discharge, Wbf is the channel width at bankfull conditions, Sbf is the channel slope at bankfull conditions and D50 is the median grain size. Using these relationships, we calculated the Froude number for bankfull conditions at a range of discharges in an idealized gravel bed river, assuming a median grain size of 10 cm and two values of the Manning's roughness coefficient. The Froude number at equilibrium conditions using the relationships proposed by [START_REF] Parker | Physical basis for quasiuniversal relations describing bankfull hydraulic geometry of single-thread gravel bed rivers[END_REF] are subcritical, typically in the range 0.2-0.3, whatever the value of the bankfull discharge. This shows that in steady uniform flow conditions, single thread gravel bed rivers are not expected to be close to supercriticality, and would not do so even after significant narrowing due to the weak dependency of Fr on channel width (eq. S3). The [START_REF] Parker | Physical basis for quasiuniversal relations describing bankfull hydraulic geometry of single-thread gravel bed rivers[END_REF] relationships (equations S4, S5) are empirical based on a dataset of 72 gravel bed rivers in the UK, Canada, and the USA, which exhibit a degree of universality between them, but do not pertain to braided rivers which are very often the state of rivers during terrace aggradation phases (e.g., [START_REF] Hanson | Fluvial terrace formation along Wyoming's Laramie Range as a response to increased late Pleistocene flood magnitudes[END_REF].

Using high resolution lidar topography and the Floodos hydrodynamic model, we explored the characteristics of the flow conditions at bankfull conditions for a braided river, the Rakaia river in New Zealand. The present day conditions of the Rakaia are a good analogue both for the initial conditions of the experiments and for the former conditions of many present day rivers that are incised with flights of terraces abandoned on the valley sides (e.g., Charwell

River, New Zealand; [START_REF] Bull | Stream-terrace genesis: Implications for soil development: Geomorphology[END_REF].

The supercritical flow conditions that initiated the development of cyclic steps in the experiments (Fig. S9) occurred after channels had narrowed following the initial stripping of the sediment cover. Therefore, we performed some simple tests using the DEM of the Rakaia river to explore the impact of different constant channel widths with a fixed discharge on the Froude number (Fig. S9A), and the impact of different discharges in a constrained fixed width channel (Fig. S9B). For a fixed discharge, a narrower channel increases the average value of the Froude number and for a fixed width, increased discharge also increases the average values of the Froude number. The average values of the Froude number remain sub-critical for the scenarios presented in Fig. S9, but the peak values of Fr approach supercriticality due to the local spatial variations in the flow conditions. The numerical simulations in Fig. S11 were carried out using a straight channel, which is overly simplistic for a natural river system. The maps of Froude number for the Rakaia (Fig. S10B) river indicate that flow conditions are not spatially uniform at bankfull discharge and are subcritical (Fr < 1) but the values of the Froude number can be higher (~0.4-0.6) than the values predicted by the relationships of [START_REF] Parker | Physical basis for quasiuniversal relations describing bankfull hydraulic geometry of single-thread gravel bed rivers[END_REF]. However, it is possible to have supercritical flow conditions during frequent floods when the flow is focussed into a narrower channel with a spatially non-uniform channel width (Fig. S10C-E).

Such conditions may be present in natural settings following a perturbation to the sediment supply or transport capacity and the subsequent stripping of sediment cover and focussing of the flow into a significantly narrower channel with spatially variable width (Fig. 2A). The simulations in Fig. S10 were carried out with no alteration to the channel slope. If the slope of the bedrock surface is steeper than the slope of the alluvial channel flowing over the

Figure 1 .

 1 Figure 1. : Example of terraces generated in the experiments. The white surface is the silica paste acting as cohesive bedrock and the brown material is sand, acting as coarse material transported as bedload. Progressive channel incision and channel narrowing leads to the abandonment of the previously active river bed, preserved as strath terraces. Inset shows topographic cross-section from A to A', with terraces labeled in yellow and orange.

Figure 2 .

 2 Figure 2. (A). Extracted long profiles during experiment 1. Green line (t = 11 min) shows the profile of the alluvial river, with the subsequent stripping of the sediment cover progressively exposing the bedrock surface within the channel. Black solid line (t = 171 min) shows the final long profile. Grey dashed lines indicate the location of preserved terraces at the end of the experiment. Instabilities form mid-channel, developing into a single knickpoint step that propagates upstream. (1): Initial alluvial channel. (2): Downstream stripping of alluvial cover. (3): Cyclic step development. (4): Knickpoint retreat in headward direction. (B): Evolution of channel width (lines) and rate of vertical elevation change (crosses) through time for the experiments (expt 1 = black, expt 2 = blue, expt 3 = red) at a cross-section 25% distance from the channel inlet. Width is plotted using a 10 point moving average and elevation change is plotted using a 5 point moving average. Dashed vertical lines indicates the timing of the perturbation in experiments 1 and 2 at 37 mins (gray), and the timing of the presence of a knickpoint (colored lines). After the perturbation, there was a reduction in channel width during both experiments 1 and 2. Focusing of the flow into a single thread occurred due to the gradual tilting in experiment 3. There was a short period of channel widening between 150 and 175 min associated with the erosion of the sediment deposited at the inlet during the initial stages, unrestricting the flow and inundating a larger area and an apparent increase in the channel width. The highest incision rates are associated with (i) the rapid downstream stripping of sediment and (ii) the headward migration of a knickpoint through bedrock.

Figure 3 .SI Section 1 :

 31 Figure 3. Conceptual model showing the stages of supercritical terrace formation. (A) Downstream stripping of sediment cover following perturbation to Qs/Qsc. (B) Under

  S1) where C if a friction coefficient in [L -1 .T -1 ], H is local water depth and S is the local channel slope. Equation S1 is obtained by integrating vertically a laminar flow velocity profile assuming steady-uniform flow. In pure laminar flow, C is entirely set by the water viscosity μ

Fig. S3 :

 S3 Fig. S3: (A) Photograph of water flowing over a flat surface (slope = 3.3 o , Q = 1 l/min), with black lines indicating edge of water extent measured manually from the image. Flow is from top to bottom in the image. (B-F) Floodos hydrodynamic model output using laminar flow option and different values of the friction coefficient. Black lines are the same as (A) and shown for comparison of modelled against observed flow extent. The predicted theoretical

Fig. S4 :

 S4 Fig. S4: (A) Photograph of experiment 1 (t = 66 mins), with wetted area indicated by the black bars. Flow is from top to bottom in the image. (B) Floodos model water mask coloured by water depth, using friction coefficient of 2.5 x 10 6 m -1 s -1 , for the same time period overlain over the photograph, with the same black bars from (A). The Floodos water mask closely matches the manual measurements of channel width from the photograph throughout the length of the channel. At the left hand side in the downstream reach, a terrace is starting to be abandoned, hence the very low flow depths in this region.

Fig. S5 :

 S5 Fig. S5: (A-C) Extracted long profiles from the three experiments. (A) Experiment 1, where the input sediment flux was stopped after 37 minutes (80 profiles in total). (B) Experiment 2, where the discharge was doubled to 3 l/min after 37 minutes (61 profiles in total). (C) Experiment 3, where the slope was increased at 1 o per 60 minutes from an initial slope of 2 o for 280 minutes (132 profiles in total). The colour of the lines indicates the evolution during the experiments. The dashed yellow lines in (C) show the initial slope of the experiment (2.5 o ), and the projected slope of the channel at the end of the tilting if there had been no erosion of the alluvial cover or bedrock (7.2 o ). The maximum slope achieved in the channel before the onset of the incision of bedrock was ~4.5 o . Grey dashed lines in (A) and (B) indicate the elevation of the bedrock surface at the start of each experiment, to show the transition from alluvial stripping and bedrock incision. The same line is not shown in (C) as the elevation of the bedrock surface changed during the course of the experiment due to the tilting. (D-F) Preserved terrace surfaces at the end of each experiment. Solid lines indicate the long profile of the active channel, and dashed lines indicate the location and elevation of terraces. The slope of the upper terraces is greater than the incised active channel.

Fig. S6 :

 S6 Fig. S6: Evolution of channel cross-sections at the 25% position from the water inlet. Crosssections were extracted from this location to ensure it was located over an area that experienced incision and the development and propagation of a knickpoint. Colours represent the same time points as Fig. S1. It can be seen that channel incision occurred under narrow channel width conditions, followed by lateral erosion of the banks leading to channel widening and minor channel aggradation in experiments 2 and 3. Grey dashed lines in (A) and (B) indicate the elevation of the bedrock surface at the start of each experiment, to show the transition from alluvial stripping and bedrock incision. The same line is not shown in (C) as the elevation of the bedrock surface changed during the course of the experiment due to the tilting.

Fig. S7 :

 S7 Fig. S7: (A) Selected long profiles showing alluvial channel just before sediment supply stopped (t = 37 min; black), channel containing cyclic steps, after sediment has been stripped away from the bedrock profile (t = 83 min; purple) and then three long profiles from different stages of knickpoint retreat (t = 102, 119, 137 min; blue, orange and red, respectively). (B) The variation in rate of elevation change with distance along the channel between the different time steps. The consistently high rate of elevation change from t = 37 minutes to t = 83 minutes (purple line) corresponds to the time period when the alluvial cover was being stripped downstream.

Fig. S8 :

 S8 Fig. S8: Exploration of flow conditions during the development of the self-formed knickpoints. (A-B) Topography during the experiments at selected time steps, with the colours indicating the water depth. (C-D) Profiles for the bed and water surfaces for the reach along the black lines in (A) and (B). (E-F) Extracted Froude number along the same reach as the (C) and (D).

  width. Substituting these two equations into each other gives the following equation for the flow depth: is calculated from the flow velocity, the acceleration due to gravity (g) and the flow depth: (Fr = U/(gH) 0.5 ). Substituting equation S2 into this equation gives the following relationship between the Froude number and the discharge, slope, width and Manning's roughness coefficient:

Fig. S9 :

 S9 Fig. S9: (A) Impact of channel width on Froude number for a fixed discharge (Q = 2000 m 3 s - 1 ). (B) Impact of discharge on Froude number for a fixed width (W = 80 m). As predicted by equation S3, both width and discharge have a small impact on the average Froude number in the channel, but the flow can become supercritical in a few locations.

ACKNOWLEDGMENTS

We thank F. Métivier for the loan of the sediment feeder, P. Davy for providing the Floodos code and T. Croissant for his help with the model runs. We thank J. Scheingross, A. Mather, A. Limaye and an anonymous reviewer for constructive comments that helped improve the final version of this manuscript. The research was funded by Marie Sklodowska-Curie Actions Individual Fellowship No. 703230 (to E.B).

Supplementary Information

sedimentary deposits (e.g. [START_REF] Finnegan | Sediment supply, base level, braiding, and bedrock river terrace formation: Arroyo Seco[END_REF], the likelihood of supercritical flow and the possible development of bed instabilities that evolve into knickpoints is increased. We therefore suggest that the processes leading to rapid strath terrace generation through selfformed knickpoint retreat, in the absence of any vertical perturbation to the bed elevation, identified in the experiments are possible in natural environments.


est destinée au dépôt

In automotive application, polymer electrolyte membrane fuel cells (PEMFCs) are operated through dynamic working conditions. Dynamic driving cycles are mainly composed by cold starting, sudden load variations, continuous load running, idle running, full power running, overload running, and start-stop repetitions [1][2]. These conditions if operated for long-time can severely impact the PEMFC lifetime. To solve this problem, a power balance strategy related to a PEMFC-battery hybrid structure is usually considered to improve the PEMFC durability [3]. As a result of the power balance strategy based on the battery state of charge optimization, a more beneficial quasi-static load profile is commonly obtained for the PEMFC operations. The aim of this work is to observe the effects of the quasi-static load profile on the PEMFC degradation. To this purpose a quasi-static reference profile was developed to separately investigate the impacts of both stationary conditions and load variations (from 0 h to 408 h) and start-stop conditions (from 408 h to 792 h) on the stack voltage degradation. A stack voltage reduction was observed after 408 h of operation, where the impact of the load variation was tested. Subsequently, the first recovering was found. Test results after 792 h also showed the start-stop impact on the stack voltage degradation.

Particularly, considering the voltage variation with respect to the cumulated energy an higher voltage slope can be observed. An accelerated load profile was finally developed to study the influence of the number of repetitions in both load and start-stop cycles on the PEMFC voltage degradation; the results are then compared. Both reference and accelerated load profiles are reported below.

Voltage

  Singularity-Based Diagnostic for H2&CO2/O2 PEM Fuel Cell Stack designed for Operation in µCHP Units Raffaele Petrone, Rungsima Yeetsorn, Fabien Harel, Daniel Hissel, Marie-Cécile Pera, Fei Gao, Elena Breaz, Stefan Giurgea To cite this version: Raffaele Petrone, Rungsima Yeetsorn, Fabien Harel, Daniel Hissel, Marie-Cécile Pera, et al.. Voltage Singularity-Based Diagnostic for H2&CO2/O2 PEM Fuel Cell Stack designed for Operation in µCHP Units. FDFC 2017, 7th International Conference on Fundamentals & Development of Fuel Cells" Jan 2017, STUTTGART, Germany. hal-01831175

FEMTO-ST, CNRS, Univ. Bourgogne Franche-Comte, rue Thierry Mieg, F-90010 Belfort Cedex, France

FCLAB, CNRS, Univ. Bourgogne Franche-Comte, rue Thierry Mieg, F-90010 Belfort Cedex, France

Department of Industrial Chemistry, King Mongkut's University of Technology North Bangkok, Bangkok, Thailand

Université de Lyon, IFSTTAR / AME / LTE, France

FEMTO-ST, CNRS, Univ. Bourgogne Franche-Comte, UTBM, rue Thierry Mieg, F-90010 Belfort Cedex, France


Philippe Meseure
email: philippe.meseure@univ-poitiers.fr

Xavier Skapin
email: xavier.skapin@univ-poitiers.fr

Emmanuelle Darles
email: emmanuelle.darles@univ-poitiers.fr

Guilhem Delaitre

Accelerated Simulation of Brittle Objects for Interactive Applications

Keywords: Physically-based animation, fractures, connected-component detection, topological models, real-time simulation

This article presents a model that aims at computing deformation and simulating fractures. To allow the use of linear elasticity in small displacements for deformation of a brittle object while still enabling any rigid motion, a rigid reference is computed using the Shape Matching method and all displacements are evaluated with respect to this reference. Fractures are handled using a stress tensor computed at each vertex of the object's 3D mesh. Some accelerations are proposed, that allow a faster determination of fracture areas and a fast processing of new connected components.

INTRODUCTION

For a long time, fracture simulation has been one of the main concerns of physically-based animation [START_REF] Terzopoulos | Modeling inelastic deformation : Viscoelasticity, plasticity, fracture[END_REF][START_REF] Norton | Animation of fracture by physical modeling[END_REF]. However, despite the amount of proposed solutions, it remains a tedious process. As far as "rigid" brittle objects are concerned, the fracture simulation consists in computing the (infinitesimal) deformation of the object when subject to external constraints and detecting where a local strain or stress exceeds some threshold. Finite element models using linear elasticity (Hooke's constitutive law) has been proposed in the pioneering work presented in [START_REF] O'brien | Graphical modeling and animation of brittle fracture[END_REF]. However, such an approach and all its derivatives require the (non-linearized) Green/St-Venant deformation tensor which makes them more computationally demanding. Real time finite element methods usually assume that displacements are small and use Cauchy's infinitesimal strain tensor [START_REF] Cotin | Real-time elastic deformations of soft tissues for surgery simulation[END_REF][START_REF] Cotin | A hybrid elastic model allowing real-time cutting, deformation and force-feedback for surgery training and simulation[END_REF]. Since, in rigid bodies, fractures occur upon small deformations, this tensor seems a good choice for the simulation of such objects. Unfortunately, it does not allow the object to rotate (it produces large deformation instead), and rigid displacements cannot directly be taken into account. This limit is crippling for rigid object simulation. In this article, we propose a model that allows both rigid motion computation and strain computation using Cauchy's linear tensor, which makes it welladapted for interactive or real time environments.

In time-critical applications, the fracture process appears as a heavy step of the simulation that should benefit from any solution that allows computation time gains. We therefore highlight a method that avoids the computation of a costly fracture criterion by ignoring places where no fracture is likely to occur. Moreover, since our method explicitly computes rigid motion, it is necessary to handle each connected component of the object separately. However, after a fracture, a simulated object can split into multiple components. Unfortunately, connected component detection is a global topological property that requires a walk through the whole mesh of the object. Ideally, this walk should be proceeded only when new connected component have actually appeared. We thus propose the use of an oracle that avoids, as much as possible, such a walk.

The paper is organized as follows: After the presentation of previous work in Section 2, the model is detailed in Section 3. Section 4 presents how fractures are handled and different process improvements. Finally, Section 5 shows the obtained results before concluding.

PREVIOUS WORK

Fracture simulation relies primarily on deformation computation [START_REF] Nealen | Physically based deformable models in computer graphics[END_REF], but many approaches exist [START_REF] Muguercia | Fracture modeling in computer graphics[END_REF]. Indeed, when a deformation is computed, different criteria can be used to make cracks appear in the structure. Geometric criteria have been proposed, such as an elongation of springs beyond a given threshold [START_REF] Norton | Animation of fracture by physical modeling[END_REF] or an excessive strain tensor value [START_REF] Terzopoulos | Modeling inelastic deformation : Viscoelasticity, plasticity, fracture[END_REF]. Force or stress thresholds can be used as well. Thus, O'Brien and Hodgins proposed a simulation of fracture based on linear elasticity and the computation of a separation tensor at each mesh vertex [START_REF] O'brien | Graphical modeling and animation of brittle fracture[END_REF]. The separation tensor of a vertex is obtained from the stress tensor of the surrounding elements. By considering positive and negative eigenvalues separately, each tensor is separated into a tensile tensor and a compressive one. Tensor components of volumes surrounding a node are associated to create a new node tensor that allows to take into account both tensile and compressive stresses. This method has been extended in various approaches, for instance [START_REF] Pfaff | Adaptive tearing and cracking of thin sheets[END_REF] to simulate 2D tearing and [START_REF] Koschier | Adaptive tetrahedral meshes for brittle fracture simulation[END_REF] for adaptive resolution models. To handle fractures, it is also possible to rely on modal analysis [START_REF] Glondu | Realtime simulation of brittle fracture using modal analysis[END_REF], but this is a global method that can be quite costly, especially for high resolution meshes. Boundary element method have also been used [START_REF] Hahn | High-resolution brittle fracture simulation with boundary elements[END_REF][START_REF] Hahn | Fast approximations for boundary element based brittle fracture simulation[END_REF]. Fractures can also be applied on meshless models [START_REF] Pauly | Meshless animation of fracturing solids[END_REF].

Cracks may be applied wherever a fracture criterion is verified. However, it is possible to focus only on high stress areas, and control cracks to avoid a shattering effect [START_REF] Koschier | Adaptive tetrahedral meshes for brittle fracture simulation[END_REF]. Local geometry is not necessarily adapted to the topological modifications that the mesh should undergo. Some methods propose to virtually cut volumes while keeping well-shaped elements but the used mesh is no longer manifold [START_REF] Molino | A virtual node algorithm for changing mesh topology during simulation[END_REF]. XFEM approaches consider discontinuous element shape functions to represent and simulate fractures [START_REF] Koschier | Robust extended finite elements for complex cutting of deformables[END_REF]. Some interactive methods tend to use predefined crack patterns that locally modify the mesh [START_REF] Müller | Real time dynamic fracture with volumetric approximate convex decompositions[END_REF][START_REF] Chen | Physicsinspired adaptive fracture refinement[END_REF][START_REF] Parker | Real-time deformation and fracture in a game environment[END_REF], but cracks are not "physicallymanaged" in such cases. The solution proposed in [START_REF] Schvartzman | Fracture animation based on high-dimensional voronoi diagrams[END_REF] appears as a compromise where the crack pattern is completely applied at the time of impact (no crack propagation) but is adapted to the deformation energy field in the fracturing body and can be user-controlled.

Different deformable models can be used to find fracture zones, but many approaches rely on continuous mechanics, solved by a finite element method, as in [START_REF] O'brien | Graphical modeling and animation of brittle fracture[END_REF] and subsequent work. Since an object can break due to small displacements, using a linearized, more computationally-efficient, strain tensor is possible, but is not compatible with rotations of the object. The corotational method consists in computing local rotations at each node to extract rotation and focus on deformation only [START_REF] Müller | Stable real-time deformations[END_REF]. This method has been used to simulate fractures (Müller et al., 2004;[START_REF] Chen | Physicsinspired adaptive fracture refinement[END_REF], but is known to produce ghost forces. To avoid such forces, a rotation can be computed, not for every node but for every element. However, in such a case, the fracture criterion should rely on the stress computed for each element (Müller and Gross, 2004) and not on each node. Indeed, a node has several possibly different displacements, one for each surrounding element. Considering a node-based measure of stress can be seen as non-consistent. Approaches based on node stress criteria for corotational FEM models have been used in some work [START_REF] Busaryev | Adaptive fracture simulation of multi-layered thin plates[END_REF] although they lack theoretical backgrounds.

To conciliate computationally-efficient deformation measures and rigid motions, some methods consider that a brittle object deforms and eventually breaks only for given situations. They propose to simulate this object as a perfectly rigid one, but switch to a static computation of (linear) deformation in case of collisions [START_REF] Müller | Real-time simulation of deformation and fracture of stiff materials[END_REF][START_REF] Bao | Fracturing rigid materials[END_REF][START_REF] Liu | Meshless simulation of brittle fracture[END_REF]. It is also the case in [START_REF] Zhu | Simulating rigid body fracture with surface meshes[END_REF] by relying only on a surface mesh to compute both deformation and crack path. Such approaches suppose that only an instantaneous impact can produce fractures. Nevertheless, to simulate a fracture resulting from a progressive deformation (see Figure 7), a dynamic approach should be preferred [START_REF] Glondu | Realtime simulation of brittle fracture using modal analysis[END_REF][START_REF] Koschier | Adaptive tetrahedral meshes for brittle fracture simulation[END_REF]. Anyway, the simulation of a rigid body simulation (with only 6 degrees of freedom) is several orders of magnitude faster to compute than a finite element mesh. In time-critical applications, the computation time of a simulation step should be kept as constant as possible. It cannot be envisaged to switch from a rigid body to a finite element model without expecting a high variation of the latency of the simulation system. Instead, for such applications, we think that simulating a breakable rigid body, at every time step and not only time step with collisions, using a low-deformable model is acceptable.

Indeed, finite element methods for small displacements appear to be fast enough to allow real-time interactions [START_REF] Cotin | A hybrid elastic model allowing real-time cutting, deformation and force-feedback for surgery training and simulation[END_REF][START_REF] Müller | Stable real-time deformations[END_REF]. To deal with rotations, rigid motions should be computed by an additional, coupled, model. The pioneering work in [START_REF] Terzopoulos | Physically based models with rigid and deformable components[END_REF] has shown how to combine rigid body laws of motion with linear deformations. Unfortunately, the equations of rigid motion and deformation are coupled and there-fore hard to solve in real time. The Shape Matching method aims at finding the minimum deformation between a given shape and a modification of this shape. This method has been used in Computer Graphics to simulate the deformation of a body by comparing the position of each node to its corresponding reference node in the undeformed body and adding a force which tends to minimize the difference [START_REF] Müller | Meshless deformations based on shape matching[END_REF]. In practice, this approach consists in inserting a spring between these two positions. The resulting simulation is fast and allows some special animation effects, but is usually not realistic. Note that, while observing a moving deformable object, the separation between the rigid motion and the deformation of the body is quite arbitrary. Shape Matching appears as a way to define, by minimizing deformation, the rigid reference component. This can appear as a nonphysical approach, since the motion of this rigid component is not guaranteed to be continuous (motion is not taken into account in the computation, only positions are relevant with this method). However, it is possible to smooth the motion of the rigid component to make it look like continuous. Shape Matching must be applied separately for each connected component of the object. More generally, any method that computes rigid motions explicitly, should precisely detect each new detected component after a fracture [START_REF] Müller | Real-time simulation of deformation and fracture of stiff materials[END_REF].

To overcome the above-mentioned limitations, this paper presents an approach that:

• Relies on linear elasticity to compute the deformation of a brittle object. Displacements are guaranteed to be small by optimizing the position of a rigid reference from which displacements are computed.

• Uses an accelerating technique based on Gerschgorin's theorem to find fracture locations rapidly.

• Uses an oracle based on "vertex splitting" to detect the appearance of any new connected component and update the computation of the rigid reference of each connected component.

AN ALMOST-RIGID DEFORMABLE MODEL

In our model, the deformation is considered to be as small as possible (remember that most brittle objects appear to be rigid at the macroscopic scale). In other words, most energy of the system should be spent for rigid motion, while the remaining energy should be used for deformation. At each time step, position and orientation corresponding to the rigid motion are both determined using the Shape Matching method and allow to define a, so-called, "rigid reference". In [START_REF] Müller | Meshless deformations based on shape matching[END_REF], the Shape Matching method was mainly used to avoid the use of a mesh, by supplying a mass/spring model with a fast way to handle its volumetric behavior. On the contrary, in our approach, the Shape Matching aims at computing the position and orientation of a rigid reference to compute a minimum deformation and maintaining, as long as possible, the hypothesis of small displacements required by the linear elasticity continuous model.

Figure 1 shows, for a block supported by immovable obstacles placed at its ends, both its deformed state and rigid reference.

Principles

Our model is a tetrahedral mesh where each vertex i is characterized by a current position x i and a rest position x 0 i (defined in an undeformed and centered configuration). At each time step, the process consists in: 1. Applying Shape Matching with the current and rest positions of the mesh vertices to obtain both position (a translation vector τ) and orientation (a rotation matrix R) of the rigid reference, 2. For each vertex, computing the reference position

x ref i = Rx 0 i + τ and its displacements as

d i = x i - x ref i , 3.
Calculating the forces applied to each vertex by each tetrahedron, using an explicit finite element approach, 4. Summing all the forces applied to vertices, inferring their accelerations and integrating them twice, 5. Handling fractures.

An explicit finite element method is used [START_REF] Müller | Real-time simulation of deformation and fracture of stiff materials[END_REF] to solve the equations of linear elasticity based on Cauchy's infinitesimal strain tensor. Let C be the 6 × 6 stress-strain matrix based on Young modulus and Poisson coefficient. For each tetrahedron (e), both 6 × 12 strain-displacement matrix B (e) and 12 × 12 rigidity matrix K (e) are computed during a pre-processing phase:

K (e) = B (e) T CB (e) (1)
At each time step and for each tetrahedron (e), the forces exerted on its vertices are computed as:

f (e) = K (e) d (e) (2)
where d (e) is a vector that gathers the displacements of the four vertices, and f (e) gathers the forces resulting from the deformation of the tetrahedron and applied on the vertices.

The use of explicit or semi-implicit integration schemes is theoretically possible, but in the case of fractures, very small time step (microseconds or tenths of microseconds) should be chosen. Indeed, to simulate fracture and avoid large deformation, a high stiffness (i.e. Young modulus) is required which directly impacts stability. Such small time steps are not compatible with interactivity time constraints. The implicit Euler method combined with the Newton-Raphson method to solve non linear systems [START_REF] Baraff | Large steps in cloth simulation[END_REF] allows large time steps but can result in a highly-damped behavior for stiff systems. Recently, exponential methods have been proposed [START_REF] Michels | A stiffly accurate integrator for elastodynamic problems[END_REF]. These approaches allow the resolution of stiff equations in a fast and robust way while not loosing as much mechanical energy as implicit methods. In our case, the implicit Euler system is solved as proposed in [START_REF] Hilde | A fast implicit integration method for solving dynamic equations of movement[END_REF]. This method is mainly dedicated to small time steps and does not require the computation of the Jacobian matrix of the system while still allowing a good control of the convergence of the system resolution. Other integration methods such as implicit or exponential ones could be used, for more efficient computation times for instance, but not with the intent of reducing the time step.

Fractures are handled as in [START_REF] Koschier | Adaptive tetrahedral meshes for brittle fracture simulation[END_REF]. Note that after one or more fractures, multiple connected components may appear. Shape Matching must therefore be applied on each connected component (R and τ are computed for each one), so a precise detection of any new component is mandatory. In our implementation, all the vertices that belong to hte same connected component share the same id. The various simulation steps that depend on connected components (for instance, the deformation computation) are slightly modified, to take into account the data related to the concerned connected component, using the id stored in vertices.

The overall algorithm for computing the behavior of a body is given in algorithm 1.

Using Shape Matching

As explained in [START_REF] Müller | Meshless deformations based on shape matching[END_REF], for a given connected component, the position of its rigid reference is defined as the position of the center of mass G of its point cloud:

x G = ∑ N i=1 m i x i ∑ N i=1 m i (3)
where N is the number of vertices in the mesh, m i the mass of each of them and x i their position. The rotation is extracted from the following matrix:

A = N ∑ i=1 p i x 0 i T (4)
where p i is the relative position of the i th vertex with respect to the center of mass, that is

p i = x i -x G .
The computation is shown in algorithm 2. Note that there is no loop such as "for each element of a connected component", which could cause computation time overheads due to dedicated walk through the topological structure of the simulated body.

This method is basically an approach that computes a position that minimizes the length of all p i (grouped in a single energy function). Unfortunately, the method ignores the history of motion. Therefore, the temporal continuity of the rigid reference's orientation cannot be ensured. In practice, when the It turns out that only isolated tetrahedra (in other words, connected components with only four vertices), exhibit such a discontinuous behavior. We have investigated a first solution that considers isolated tetrahedra as rigid, since they are atomic elements of the initial mesh. This can be done in our algorithm using only a slight modification, by considering that rigid bodies can be simulated using a constrained deformable model [START_REF] Van Overveld | All you need is force: a constraint-based approach for rigid body dynamics in computer animation[END_REF]. More precisely, the vertices are first moved according to external forces, then Shape Matching is computed for the tetrahedron, as done for any connected component of the model. After that step, current vertex positions x i in isolated tetrahedra are projected onto their corresponding reference position (x ref i ). Velocity is adapted accordingly. This solution has a drawback: Isolated pseudo-rigid tetrahedra tend to stop rotating quite quickly. Indeed, angular momenta are not explicitly maintained by this approach, rotations result from the difference of velocities of the vertices of the tetrahedra. Since vertices are forced to keep their tetrahedron undeformed, they tend to loose velocity and therefore, the orientation of their tetrahedron steadies rapidly.

We have investigated another solution to compute rotation of isolated tetrahedra, by using a QR decomposition as proposed in [START_REF] Nesme | Efficient, physically plausible finite elements[END_REF]. This method provides better stability for tetrahedra (i.e. rotation matrix continuity), provided that the QR decomposition is handled with the same vertex order between two steps.

Display

Even if a high Young modulus is applied, the object can still exhibit visible deformations due to the computation time step which is too large for fractures (see Figure 1). To guarantee an accurate display of the object, the rigid reference can be displayed instead of the deformable body itself. This approach works well in practice, but has an impact on the collision process of the simulation, if overlaps must be avoided. Fortunately, since the reference body is not physically simulated, any modification of position and orientation can be applied arbitrarily, for display only. Any correction method that controls the position and orientation of a rigid body, such as an optimization algorithm to minimize overlap, can therefore be applied.

In our simulation, no correction was however necessary, since overlaps were low enough in practice. Shape Matching is applied independently on each connected component. Should a new connected component appear after a fracture, it is unlikely that its newly-computed position and orientation are continuous with respect to the position and orientation it had previously, when it was still connected. Figure 2 illustrates this discontinuity problem. To avoid the discontinuity of display, an interpolation between its previous state and its newly-computed one is applied. More precisely, a current state is first computed, based on its original connected component. This can be easily made by applying Shape Matching on x ref i positions. Next, a target state is computed by applying Shape Matching on {x i } positions. After a fracture, current and target states differ. Target state should be used for the display of its corresponding connected component, but its position may be too far from the position where it was previously located. Current state corresponds to the location of the connected component at the fracture time. During the following time steps, we choose to modify the position and orientation of the current state so that it tends to reach the target state. Although a proportionalderivative controller could be used for this purpose, only a simple proportional controller has been used in our implementation and gives satisfying results in practice. Let τ cur and τ target be the current and the target translations respectively. Then τ cur is updated to τ cur + α(τ targetτ cur ), where α is a small value (typically 0.1 or less). A similar approach is used for rotation, based on quaternions, since they offer a convenient way to interpolate rotations [START_REF] Shoemake | Animating rotation with quaternion curves[END_REF]. If q is a quaternion representing the rotation that transforms the current orientation into the target orientation of a connected component, the quaternion q + 1 represents (after renormalization) half this rotation. By repeating this simple formula several times, we can quickly divide the angle of rotation between current and target orientation by any power of 2. Accordingly, the α value used for translation is chosen as the corresponding negative power of 2.

Note that the target state is also updated at each time step, since all x i positions change. A more general approach, where a rigid reference is simulated as a perfect rigid body, linked by two springs, one linear and the other one angular, to the deformed configuration, and subject to collisions constraints is another way to compute the current state at each time step. This approach is physically-based but results in a 2 ndorder control of the current state used for display, so possibly an undesirable latency.

FRACTURE PROCESS

Due to the linear interpolation functions used by the finite element method, the stress tensor is constant all over each tetrahedron and can be computed as:

σ (e) = CB (e) d (e)
(5) As proposed in [START_REF] Koschier | Adaptive tetrahedral meshes for brittle fracture simulation[END_REF], a node stress tensor can be computed for each vertex as:

σ = ∑ (e) m (e) σ (e) ∑ (e) m (e) (6)
where m (e) is the mass of the element (e). Here, the expression has been slightly modified compared to the original article, but the result is the same. Our notation only aims at showing that it is basically a weighted average value of surrounding element stress tensors. Note that, if each tetrahedron provides each of its vertices with a quarter of its mass (mass lumping technique), the denominator of equation 6 corresponds to the mass of the concerned node multiplied by 4. A fracture is applied on some vertex when at least one eigenvalue of its stress tensor exceeds a given threshold ξ. A fracture plane is built using the related eigenvector as its normal. Each surrounding tetrahedron is marked as "positive" or "negative" depending on the side of the plane where its center is positioned. Any face that is adjacent to a pair of tetrahedra with different signs is split to separate them. The whole process aims at splitting a fan of faces around the "fracture vertex" to make it split. If a fracture generates new connected components, these components must be identified to allow distinct Shape Matching applications when needed. This represents the last step of the fracture process.

To accelerate this process, two lines of approach have been considered: a fast location of fracture zones and a fast identification of new connected components. The algorithm 3 give an overview of the fracture process.

Finding fracture location

One important time loss during fracture location comes from the necessity to compute eigenvalues of stress tensors, merely to determine if one of them may exceed some threshold ξ. In mechanical engineering, this step is often sped up using the Gerschgorin's theorem which, to our knowledge, has never been proposed to the Computer Graphics community. This theorem states that an upper bound of the eigenvalues of a square matrix M can be computed considering the sum of absolute values of some of its terms. In our case, M is a symmetric matrix, so the upper bound can be found by computing the maximum value of three terms b i , one for each column i:

b i = M ii + ∑ j =i M ji (7)
After calculating the stress tensor of a vertex, the associated b i values are computed. If all of them remain below the fracture threshold, no fracture appears at this vertex and the eigenvalues computation can be avoided. On the contrary, if the test fails, the eigenvalues must be computed since, if at least one of them is greater than the given threshold, its associated eigenvector must be known for the following steps. Note that the overall complexity is small, since only 6 floating point additions, 6 absolute values and 3 comparisons are necessary to complete Gerschgorin's test at a vertex.

Connected Component Determination

Each connected component of the body is characterized by a rigid reference, defined as a current state (position and orientation) and a target state. When a fracture occurs, one or more new connected components might appear, which should be clearly identified to update their states. Connected components identification requires to scan all the elements (for instance vertices) of the current connected component and mark them (more precisely, in our implementation, "marking vertices" means providing them with the id of their connected component). If some elements are not marked (i.e. these elements still store the id of their connected component before fracture), that means that another connected component must be identified. Reversely, if all elements are marked, that means that no other connected component exists. This is a time-consuming process (linear with respect to the number of elements to be marked, i.e. vertices or volumes) that should be avoided when no new connected component is created (in other words, when a first walk would uselessly mark every element).

We propose to use an oracle that can test this case, in order to speed up the overall component identification. Of course, it must be efficient enough to avoid false positive cases as much as possible (that is, cases where the oracle announces new connected components whereas no such component has appeared actually). Moreover, the oracle should be fast enough: If its computation algorithm is too complex, it will not be able to compete with a complete walk through the structure. We found an efficient oracle by considering how many vertices split after each face split. Indeed, when two adjacent tetrahedra are separated, the support edges and vertices of the face must necessarily split to allow the creation of a new connected component. In Figure 3, it is shown that the three vertices and three edges of a triangular face actually split when a new connected component appears. Note that if three vertices split, then the three edges necessarily split. Therefore, after a face split, the vertices and only the vertices of the face are checked for a split. If, at least, one of these does not split, it is ensured that no new connected component has been created. On the contrary, if the three vertices split, then a new connected component may appear, and this should be checked further. Note that the vertex split criterion must be checked after the split of each face of a fan and not once the overall split is over. An appropriate topological model should be used to ensure that the vertex split is tested effectively [START_REF] Lienhardt | N-dimensional generalized combinatorial maps and cellular quasimanifolds[END_REF]. Note that both criteria are local whereas the connected property is global. It cannot be expected that our oracle produces only right predictions, false positives are still possible: For instance, such a phenomenon appears when two parts of an object are locally separated, but are still joined if the entire structure is considered, as shown in Figure 4. When using the vertex split oracle, face ordering may sometimes influence the final result and produce false positive cases,as shown in Figure 5. For each subfigure, the red arrow shows the last face that is split. In (a), only two vertices (green and yellow) split and in (b), only one vertex (yellow) splits. In such cases, the oracle is right and states that there cannot be a new connected component (which is the case: the surrounding parts remain adjacent, for instance around the black vertices, since these vertices do not split). In (c), how-ever, the last face makes three vertices split (red, green and yellow ones). In this case, the oracle is false, since it concludes that a new connected component has appeared. Note that, if only one other face remained to be split at that time, the central vertex of the fan would not split (it splits only when the last face is split) and the oracle would not conclude that a new connected component appears. As a consequence, faces of fans with an edge at the surface of the object (as it is the case in the figure) should be checked first if possible, because these faces produces false positive if split at the end of the process. The algorithm 4 show the complete split process.

RESULTS

Several examples have been simulated and are represented in Figures 6 and7. In both cases, the simulated object is the same. The first one (Figure 6) corresponds to a classical fracture scenario where rigid objects are shot toward a breakable one. This example shows the shattering effect obtained when a fracture is applied wherever the fracture criterion is true. The second case study (Figure 7) corresponds to a tension that progressively deforms an object so that it eventually breaks. First, this example shows that fractures may appear well after a collision has occurred, so the dynamic approach of fracture is mandatory. Second, it illustrates that the display of the rigid reference is very useful. Indeed, in practice, the deformable mesh undergoes visible deformations before fracturing (this is the reason why the object does not break immediately). The display of the deformable mesh would not have produced the same visual effect. Another case study aims at simulating a more complex object, namely the Stanford Bunny (Figure 8). The mesh parameters used by the simulation are given in table 1. c), if the pointed face is the last face to be split, at least one of its vertex (in black) does not split. In (d), if the pointed face is the last face to be split, the oracle predicts a new connected component as its three vertices (red, yellow and green) split. This face, with an edge at the surface of the mesh, should be split when at least one other face of the fan has not been split yet.

In these scenarios, collision and self-collision detection has been handled using an approximative but fast approach. They were mainly based on a penalty method: Any environment obstacle is approximated by a force field (for instance, ground or rigid bodies). Every tetrahedron of the breakable body is bounded by a sphere that also produces a repulsive force field for points located inside it. Each vertex is checked against all the force fields (a regular grid is used for speeding up this step). Collision between vertices and tetrahedron force fields are taken into account only if they do not belong to the same connected component (this is quickly determined as each vertices knows the id of its connected component). A more precise but efficient collision detection and response could be used [START_REF] Glondu | Fast collision detection for fracturing rigid bodies[END_REF], but requires to be adapted to deformable bodies if used in our approach.

Note that the Broyden method to solve the implicit integration equations does not converge as quickly as a Newton-Raphson resolution. Note also that no parallelization (CPU multithreading or GPU) has been applied on the algorithms, except the use of SIMD processor instructions for FEM computations. We are completely aware that these times could be largely improved as most similar approaches are usually based on GPU programming. However, the frac- ture process is hard to parallelize in practice (because of conflict access to the topological structures of the object). This is why we did not worry on parallelizing the simulation itself to get a better performance and chose to focus on fracture only.

The average and maximum simulation times are shown in Table 2 for both corotational and our shapematching-based method. These times correspond to a complete time step, including collision and self-collision detection. The same object with the same physical properties are simulated in both methods. As shown, the corotational method is always slower than ours. This mainly shows that the rotations for each volume element always require more time to compute than a shape-matching on the connected components. Applying fractures on a simulated model usually implies a computation time overhead, about 20% for simple models but that can be up to 50% for the stanford bunny (which involves a lot of multiple synchroneous fractures). The maximum computation time when applying fractures is also mentioned in the table, although we do find this time very significant. Indeed, we have noticed that this maximum time usually occurs once, at the beginning of a fracturing phase, and is not reached anymore during the subsequent fractures.

Table 2: Average and maximum simulation times obtained for our different scenarios using corotational and our method without fracture (i.e. with an unreachable fracture threshold). For our method, the average and maximum computation times with fractures is also given. All these times include collision and self-collision detection. The effect of the different accelerations have been evaluated. First, the impact of Gerschgorin's test has been measured. This test is, in average, 3 times faster than an immediate computation of the eigenvalues. Table 3 shows the effectiveness of this test, since about 80% or more successful tests actually correspond to a large eigenvalue that implies a fracture.

The split vertex oracle for new connected components has also been tested. In table 3, the efficiency, that is the ratio of true positives, for the split-vertex oracle has been measured, in our different test cases. It can be seen that, for over 90% of cases, the oracle is right. So the oracle is a satisfactory way to detect if some new component has appeared. Note also that, in all our experiments, the oracle has never produced any false negative.

Table 4 shows the computation times dedicated to the fracture process. Several comments can be made. First, the worst case for the basic rupture process, that, by default, searches for new connected components, occurs when no new connected component is created. Indeed, all the structure must be (uselessly) walked through to find this result. This proves that an oracle that speeds up the new connected component detection is really of interest. Next, the additional cost of using split-vertex oracle is measured by comparing the computations times in case of new connected component creation using the oracle or not. This cost represents of at most 2% of computation time. But, this cost is clearly justified, because, when no new component appears, the rupture time can be reduced by, at least, 55%, and the global reduction can reach 67%. To provide a meaningful comparison, we have implemented the same node tensor fracture criterion on a corotational FEM (even if, in previous work section, we pointed out that this approach was lacking a theoretical background). This implementation also uses the Gerschgorin's theorem to speed up the broadphase that seeks for fracture nodes. We found that corotational FEM tends to produce lower node stress tensors (which is quite expected: Node displacements are lower since element rotations are cancelled). As a consequence, the shattering effect is reduced and the final animation is quite different from the one obtained with our method based on Shape Matching.

LIMITATIONS

The obtained fracture results heavily depend on the initial tetrahedral mesh of the simulated object. Remeshing techniques could be used to adapt the mesh in the fracture zone, as several approaches have already proposed [START_REF] Su | Energy stability and fracture for frame rate rigid body simulations[END_REF][START_REF] Müller | Real time dynamic fracture with volumetric approximate convex decompositions[END_REF]. These techniques are compatible with our method, provided that the rest positions and mass of the new vertices that appear after a remeshing are determined and the mechanical data of elements (stiffness and strain/deformation matrices) are updated accordingly, which can be done at the price of a little computation time. These techniques apply on a given connected component, so no side-effect on existing connected components are expected during the remeshing.

We actually investigated a local remeshing solution based on tetrahedron subdivisions based on Loop scheme as proposed in [START_REF] Meseure | Adaptive resolution for topology modifications in physically-based animation[END_REF]. This kind of remeshing is adapted to physical simulation since the shape of most elements is preserved, although scaled by a factor 1/2. This solution works well in practice for deformable body simulation, however, it is not adapted to fracture. Indeed, after a subdivision, the strain-displacement matrix B (e) of each new elements is scaled by a factor 2 and its volume, and consequently its mass, by a factor 1/8. The node stress tensor (Equation 6) is directly impacted by such modifications and the choice of mass as a weighting coefficient is no longer adapted if the resolution of the mesh is not uniform. Without a meaningful computation of node stress tensors, the criterion used to detect fracture (a high eigenvalue) is no longer relevant. As proposed in [START_REF] Pfaff | Adaptive tearing and cracking of thin sheets[END_REF] in 2D, a more appropriate tensor, that is, independent on mesh resolution, should be computed. More investigations are needed.

CONCLUSION AND FUTURE WORK

This paper has presented a model that allows the computation of small deformations of a brittle object while still allowing rigid motions. To speed up fracture management, the Gerschgorin's test has been proposed to easily locate fracture zones. A simple oracle is also used to prevent the execution of heavy algorithms needed in case of new connected components. These improvements provide interesting computation times, although no parallelization has been exploited to compute deformations. It is definitely one of our major goals, although the fracture process itself is not really parallelizable. We also want to improve visualization, by displaying more complex fracture surfaces as done in [START_REF] Koschier | Adaptive tetrahedral meshes for brittle fracture simulation[END_REF] (it would only be a visual artifact, with no cost on the fracture process itself). Furthermore, one remaining problem concerns the number of cracks (that is fan splitting) that can be handled in one step. Indeed,

Figure 1 :

 1 Figure 1: Model/reference coupling: A deformable model (wireframe) is computed using dynamics laws and a rigid reference (filled faces) is computed using Shape Matching.

Figure 2 :

 2 Figure 2: Discontinuity issue with Shape Matching. (a) A red block is deformed (red dotted line) and Shape Matching provides a position of the rigid reference (solid line). (b) A fracture occurs and the object is split. Shape Matching is separately applied on both connected components and results in a discontinuous state with respect to their positions in (a).

Algorithm 3 :

 3 findFractures() // Compute each tetrahedron stress tensor and report it to its vertices foreach Tetrahedron t do Vec12 d = getVertexDisplacements(t) ; t.σ = matC × t.matB × d ; foreach Vertex v of t do v.σ+ = t.mass × t.σ ; end // Average each vertex stress tensor foreach Vertex v do v.σ/ = (4 × v.mass) ; end // Check where a fracture occurs foreach Vertex v do if gerschgorinTest(v.σ, ξ) then Vec3 eigenval = eigenValues(v.σ) ; if max(eigenval) ≥ ξ then split(v, eigenvec(max(eigenval)) ; end end

Figure 3 :

 3 Figure 3: Oracle to detect if one new connected component appears after a face split. In (a), only one vertex and two edges split, so there is no new connected component. In (b), three vertices and three edges split, so one new connected component may appear.

Figure 4 :

 4 Figure 4: False positive of any local oracle. A local separation has occured, so any local oracle predicts a new connected component, but the separated faces still belong to the same connected component (initially topologicallyequivalent to a torus).

Figure 5 :

 5 Figure 5: False positive of split vertex oracle. (a) represents a series of adjacent edges (in magenta) at the mesh surface (in blue). A fan (dotted lines) is defined around a vertex (in yellow), inside the mesh, just below the surface faces. In (b) and (c), if the pointed face is the last face to be split, at least one of its vertex (in black) does not split. In (d), if the pointed face is the last face to be split, the oracle predicts a new connected component as its three vertices (red, yellow and green) split. This face, with an edge at the surface of the mesh, should be split when at least one other face of the fan has not been split yet.

Figure 6 :

 6 Figure 6: Two balls striking a monolith.

Figure 7 :

 7 Figure 7: A brittle block collapses under its own weight.

Table 1 :

 1 Different test scenarios and their characteristics

	Scenarios	#Tetras #Vertices
	Struck monolith	413	174
	Strained block	413	174
	Falling bunny	6142	1363

Table 4 :

 4 Computation times for the rupture process in different scenarios with and without split-vertex oracle, for ruptures that create new connected components (CC) or not.

		Without Oracle (ms)	With Oracle (ms)
	Scenarios	new CC	no new CC	new CC	no new CC
	Struck monolith	0.97	1.2	0.99	0.5
	Strained block	0.96	1.06	0.98	0.48
	Falling bunny	30.73	32.3	31.2	10.66

Table 3 :

 3 Efficiency of Gerschgorin's test and split-vertex oracle for different scenarios. The third column represents the number of fan splits. The fourth column represents the number of fractures that actually result in the creation of new connected components (CC) and the corresponding ratio. The fifth column represents the number of cases where the split-vertex oracle predicts new connected components, the corresponding ratio, as well as the efficiency calculated by checking if new connected components actually appeared. crack is allowed per simulation step, no shattering effect occurs in general. The number of handled cracks should however be bounded.

		Gerschgorin #fractures #CC creations	#new CC predictions
	Scenarios	efficiency		#	ratio	#	ratio	efficiency
	Struck monolith	85%	198	112 56.6% 119 61.1%	94.1%
	Strained block	92%	190	104 54.7% 111 58.4%	93.7%
	Falling bunny	78%	1639	759 46.3% 833 50.8%	91.1%
	if only one

ACKNOWLEDGMENTS

We would like to thank the anonymous reviewers for their comments that helped us to improve this paper. This work has benefited from the financial support of a M.Sc internship from the MIRES federation. Figure 8: Shattering fracture of Stanford Bunny. The bottom picture aims at highlighting the fracture surface and the different connected components.