Corinne Lampin-Maillet
email: corinne.lampin@cemagref.fr

Marielle Jappiot

Marlène Long

Denis Morge

Jean-Paul Ferrier

Characterization and mapping of dwelling types for forest fire prevention

Keywords: Spatial analysis, dwelling typology, GIS, forest fire risk, prevention, wildlandurban interface

, and develops a method based on GISgeo treatments to characterize different types of dwelling with regard to fire risk. Three types of dwellings were classified: isolated dwellings, scattered dwellings and clustered dwellings, using criteria based on the distance between houses, the size of clusters of houses and housing density, which can be mapped automatically. Within dwelling types, the density value of forest fire ignition changed and was twice as high for isolated dwellings as for clustered dwellings. The spatial organization of dwellings seems to have a real impact on fire occurrence. Thus maps of different dwelling types can be interpreted for use in developing fire fighting strategies or prevention actions concerning end-users such as forest and land planning managers or firefighters.

1.

Introduction

Forest fires affect extensive areas and cause serious damage that has ecological, social and economic consequences leading to very high costs of prevention and fire-fighting. Every year more than 2,600 forest fires affect about 24,200 hectares of vegetation in the French Mediterranean region (Promethée database, 1973(Promethée database, -2007)). Many large fires are linked with the dramatic land transformation that has been taking place in the Mediterranean region for some decades and that is increasing the risk of forest fire. On the one hand, agricultural fallows and orchards are slowly colonized by vegetation, while on the other hand, forest is not sufficiently utilized, where both result in the increased accumulation of fuel load. In addition, especially in the South of France, urbanization coupled with the phenomenon of forest extension is generating new spatial configurations called wildland-urban interfaces (WUI), where WUIs are commonly defined as areas where urban areas meet and interact with rural lands [START_REF] Vince | Forests at the wildland-urban interface: conservation and management[END_REF] creating a new conjunction of housing and vegetation characteristics [START_REF] Stewart | Defining the Wildland-Urban Interface[END_REF].

In recent decades, there has been a significant increase in housing all around the world [START_REF] Stewart | Defining the Wildland-Urban Interface[END_REF]Radeloff et al., 2005a). Urban sprawl observed in or near forests represents a significant threat to the environment [START_REF] Johnson | Environmental impacts of urban sprawl: a survey of the literature and proposed research agenda[END_REF][START_REF] Radeloff | The wildland-urban interface in the United States[END_REF]. Indeed, it leads to habitat fragmentation [START_REF] Theobald | Estimating the cumulative effects of development on wildlife habitat[END_REF], threatens wildlife populations by introducing exotic species [START_REF] Soulé | Land-use planning and wildlife maintenance: guidelines for conserving wildlife in an urban landscape[END_REF], and results in a decrease in biodiversity [START_REF] Mckinney | Urbanization, biodiversity, and conservation[END_REF]. Alavalapati et al. (2005) summarized the different socioeconomic and environmental concerns in a special issue. However, the main concern in WUIs is wildland fire [START_REF] Davis | The wildland-urban interface: Paradise or battleground?[END_REF]. According to [START_REF] Stewart | Defining the Wildland-Urban Interface[END_REF] "the term wildland-urban-interface is now used almost exclusively in the context of wildland fire". The WUI is becoming a priority area for fire prevention and suppression [START_REF] Stephens | Forest fire causes and extent on United State Forest Service lands[END_REF]. Indeed, there is an increasing concern about fuel reduction [START_REF] Conard | Changing fuel management strategies -The challenge of meeting new information and analysis needs[END_REF], removal of fuel in the immediate vicinity of homes [START_REF] Cohen | Preventing disaster: Home ignitability in the wildland-urban interface[END_REF] and plans for emergency evacuation [START_REF] Church | Mapping evacuation risk on transportation networks using a spatial optimization model[END_REF]. In this way, WUIs have become the central focus of wildland fire policy in Europe as well as in the United States.

We hypothesize that fire behaviour in WUIs is influenced by the pattern of urban areas within a natural landscape. Thus, it is first necessary to qualify and discriminate the spatial pattern of housing units in order to define different classes of houses: isolated dwellings, scattered dwellings or clustered dwellings. There is little literature in Europe about characterization of dwelling types through quantitative criteria in regards to fire risk assessment. For the purposes of this paper, we define 'isolated dwellings' as areas that are not very urbanized and located in a natural space. 'Scattered dwellings' refers to houses spread out over an extended urban area, and 'clustered dwellings' refers to areas where the population is concentrated. Taking into account background methods developed to characterize towns, we develop a method to characterize and stratify urban areas using criteria such as distance between houses, housing density, etc. After a brief synopsis of results stemming from an exploratory process for the characterization of dwelling types, the main aim of the paper is to present a method with which to characterize different types of dwelling with regard to fire risk. This research is a step towards defining the relationship between specific spatial organisation in WUIs and fire behaviour. We use the most efficient method to determine the typology of dwellings and combine this information with the spatial distribution of forest fire ignition points.

Research issues on urban area characterization

Many studies use land cover data, which can accurately capture high density residential areas but does not work very well in low density residential areas [START_REF] Ridd | Exploring a V-I-S (vegetation-impervious surface soil) model for urban ecosystem analysis through remote sensing : Comparative anatomy for cities[END_REF]. Census data have been widely used to measure urban sprawl in the United States [START_REF] Theobald | Land-use dynamics beyond the American urban fringe[END_REF]. Census data are collected during surveys and large scale data are easily accessible. According to [START_REF] Stewart | Defining the Wildland-Urban Interface[END_REF], housing density is the most appropriate metric for human presence. In census block groups, housing density has been used to identify WUIs across the US [START_REF] Stewart | Characteristics and location of the Wildland-Urban-Interface in the United States[END_REF]. However, there is still a need for thresholds in defining urban, WUI, and wildland areas on the basis of housing density. Depending on the definition provided by different authors, different classes of housing density are reported [START_REF] Theobald | Land-use dynamics beyond the American urban fringe[END_REF][START_REF] Stewart | Characteristics and location of the Wildland-Urban-Interface in the United States[END_REF][START_REF] Stewart | Defining the Wildland-Urban Interface[END_REF][START_REF] Hammer | Characterizing dynamic spatial and temporal residential density patterns from 1940-1990 across the North Central United States[END_REF][START_REF] Curtis | Characteristics of the private nuisance wildlife control industry in New York[END_REF]. [START_REF] Kamp | Using GIS to identify potential Wildland-Urban Interface areas based on population density[END_REF] used population density instead of housing density to define the potential WUI area in the middle of a density spectrum, which varies from 40 to 400 people per square mile. However, according to [START_REF] Alberti | the impact of urban patterns on aquatic ecosystems: an empirical analysis in Puget lowland sub-basins[END_REF], these metrics do not discriminate well between different landscape patterns, essentially because the number of housing units has increased at a faster rate than that of the human population [START_REF] Lepczyk | Spatiotemporal dynamics of housing growth hotspot in the North Central U.S from 1940 to[END_REF][START_REF] Stewart | Defining the Wildland-Urban Interface[END_REF]. [START_REF] Theobald | Expansion of the US wildland-urban interface[END_REF] developed a method for mapping the WUI at the national scale using housing density but indicated that for planning and management activities, more detailed data are necessary. Buildings should be located based on aerial photographs or using field-based data, but this kind of data is not available in all countries. Various methods corresponding to the generalization of spatial data have been developed to characterize towns: the envelop convex method, Delaunay triangulation, and recently a generalization pattern from the BD TOPO® from the French National Geographical Institute database [START_REF] Boffer | Méthode de création d'information multi niveaux pour la généralisation cartographique de l'urbain[END_REF][START_REF] Ruas | Modèle de généralisation de données géographiques à base de contraintes et d'autonomie[END_REF]. The latter method characterized towns as areas that contain similar numbers of buildings, which are merged before simplifying their shape using the mathematical morphology approach developed by [START_REF] Matheron | The birth of mathematical morphology[END_REF] and also considers recent studies from [START_REF] Edwards | Preserving the pattern of density in urban network simplification[END_REF]. Other authors [START_REF] Le Corre | Intérêts et limites de l'utilisation des images satellites en vue de la constitution d'un SIG pour la région urbaine d'Istanbul[END_REF][START_REF] Kalhori | Sectorisation de la qualité de l'air[END_REF] use mathematical morphology to delineate towns or data aggregation methods based on length criteria.

All these methods are based on specific distances. Some distance concepts are used to discriminate types of dwellings. However, there are no clear guidelines on what distance to use. A buffer distance should include a range of different fire-fighting objectives. [START_REF] Theobald | Expansion of the US wildland-urban interface[END_REF] assumed an 80 metres radius "firewise" treatment zone around each housing unit (defensive space). [START_REF] Cohen | Preventing disaster: Home ignitability in the wildland-urban interface[END_REF] proposed a distance of 40 metres as "the likelihood of it being destroyed in a wildfire". In understanding which factors influence building density and fragmentation, [START_REF] Gonzalez-Abraham | Building patterns and landscape fragmentation in northern Wisconsin, USA[END_REF] applied 50, 100 and 250 metres disturbance zones around each building. [START_REF] Stewart | Defining the Wildland-Urban Interface[END_REF] chose a 1.5 mile buffer distance based on the estimation of how far a firebrand can fly. Other methods combine both density and spatial patterns of buildings [START_REF] Theobald | Estimating the cumulative effects of development on wildlife habitat[END_REF]. Despite the significance of the WUI concerns, particularly in land and fire management [START_REF] Radeloff | The wildland-urban interface in the United States[END_REF][START_REF] Zhang | The importance of scale in using hierachical Census data to identify the Wildland-Urban interface[END_REF], empirical data on its extent and location are lacking.

An exploratory process to identify types of dwelling

The first approach was to assume that dwelling types can be discriminated according to housing density values. One easy method was to use the existing density function of the Spatial Analyst extension of Arcgis® 9.1 software [START_REF] Lampin | Dwellings characterization and mapping for forest fire risk prevention[END_REF]. The size of the radius corresponds to the extension of the wildland-urban interfaces (Lampin et al., 2006a[START_REF] Lampin | WUI and road networks/vegetation interfaces characterizing and mapping for forest fire risk assessment[END_REF].The result is relatively easy to obtain because of the existing function in the software, which is an advantage, but the cartographic figuration is not completely satisfactory. Indeed a house and its surrounding 100 metres (wildland-urban interface delimitation) was not always classified as belonging to only one type of dwelling, even when we did not use kernel density. In fact, kernel density generally increased the problem of this phenomenon. Part of the area concerned may belong to an isolated dwelling class, whereas another part of the same area may belong to the scattered dwelling class. In addition, calculation of housing density is obviously intuitive. The difficulty lies in particular, in the need to determine the best density thresholds that would allow us to distinguish the different types of dwellings. This depends on urban management decisions that may differ with the municipality. This is why we suggest that threshold values should be determined through analysis of the spatial distribution of housing.

The second approach was to consider specific distances between houses to create density classes based on the spatial distribution of houses and to discriminate types of dwellings. Rules referring to basic length concepts relating houses to the urban continuity were set, where these rules were based on existing distance concepts. The French National Institute for Statistics and Economic Studies (INSEE) defined the concept of "urban unit" as all houses located less than 200 metres apart. A distance of 200 metres can easily be covered on foot by anybody who wants to go from one house to another. Beyond 200 metres, a house is considered as an isolated dwelling. The French National Geographical Institute (IGN) created the "urban area" concept [START_REF] Cete Lyon | Les bases de données géographiques d'occupation du sol : Volet tâche urbaine[END_REF] in its geographical database BD CARTO®. These urban areas are restricted to aggregated areas which correspond to houses located less than 100 metres apart. Using this "urban area" concept, we defined clustered dwellings as houses located less than 100 metres apart. In addition, the generalization pattern from the BD TOPO® (IGN) database [START_REF] Boffer | Méthode de création d'information multi niveaux pour la généralisation cartographique de l'urbain[END_REF][START_REF] Ruas | Modèle de généralisation de données géographiques à base de contraintes et d'autonomie[END_REF] identifies the centre of town by clusters of houses that are less than 50 metres apart. The number of houses can also be taken into account; urban rules consider that houses are clustered when a group contains at least two houses separated by less than 100 metres. The method based on mathematical morphology principles was performed using IMAGE J software. Using specific distances between houses to identify the types of dwellings, the basic mathematical morphology operations (closing and opening) enabled us to discriminate three classes of dwellings [START_REF] Lampin | Dwellings characterization and mapping for forest fire risk prevention[END_REF]. This approach provides mean values of housing density based on the spatial distribution of houses which can be used to create thresholds for existing methods or to characterize each type of dwelling. But small clusters of 2 or 3 houses can still appear as isolated dwellings despite a short distance between houses. Regarding the risk of fire, this configuration of houses is as vulnerable as an isolated house. Thus the classification of types of dwellings based only on the house density criterion based on distance between houses is not adequate and requires another approach.

4.

Materials and method

Study area

The study area (fig. 1) is located in south-eastern France near Aix-en-Provence and comprises 64,700 hectares corresponding to 28 municipalities. In this area, wildland-urban interfaces are very common and urban sprawl is replacing agricultural fallows near or within forested areas. The risk of forest fire is particularly high.

Data, software

This study used data from the French National Geographical Institute (IGN), housing data from the BD TOPO® database produced in 1992, and aerial photos from the BD ORTHO® database produced in 2003, which have a resolution of 0.50 metres. Concerning housing data, the vector layer of houses was extracted from the 1992 BD TOPO® and updated from the 2003 aerial photos derived from the BD ORTHO®. The study is only concerned with residential buildings (dwellings) and includes 50,200 houses on the study site. A digitalized database of fire ignition points was used. This database has been created by the National Forest Institute and counts fire ignition points corresponding to fires that spread over more than one hectare during the period 1996-2007. 227 ignition points are contained in the study area.

One computer program was used to characterize housing types. Arcgis® 9.1, developed by ESRI (Environmental Systems Research Institute), enables the user to visualise, explore, and analyse spatial data. It is suitable for vector or raster datasets. We mainly used it to carry out geo-treatment, to combine different layers and to map results.

Identification of dwelling types using distance and counting criteria

The approach was to consider specific distances between houses combined with housecounting criteria. Rules determining basic length concepts between houses and numbers of houses were defined based on existing distance practice in the context of fire prevention. A method based on classical geo-treatment usually performed in Geographical Information System (GIS) was developed. Isolated dwellings (ID) refer to houses (or clusters of 2 or 3 houses) located more than 100 metres apart. This number of houses takes into account fire prevention management; it is difficult to protect houses against fire without endangering fire fighters. Scattered dwellings (SD) correspond to clusters of 4 to 50 houses located less than 100 metres apart. This number of houses refers to a fire-fighting strategy. A fire-fighter team composed of four trucks can protect up to 10 houses if they are located relatively close to one another, less than 100 metres apart. Clustered dwellings (CD) refer to clusters of more than 50 houses located less than 100 metres apart. In this last category, two sub-classes have been determined: a first one identifies dense clustered dwellings (DCD) corresponding to clusters of less than 10 houses located less than 30 metres apart and a second one identifies very dense clustered dwellings (VDCD) corresponding to clusters of more than 10 houses located less than 30 metres apart that being equivalent to well urbanized areas.

Different operations were carried out with Arcgis® 9.1 software. Geo-treatment enabled us to create buffers with different diameters (50 metres and 15 metres). The overlaying buffers were merged and the houses within each buffer zone were counted to differentiate housing types according to the rules presented below. For each house of a different class (isolated, scattered or clustered) a buffer of 50 metres was created. To avoid overlaying of buffers, a rule of predominance was established, where the buffer corresponding to the most clustered class took priority.

Calculation of spatial criteria related to the risk of forest fire and statistical analysis

Spatial criteria related to the risk of forest fire were determined for each type of dwelling. Calculations were made for each polygon created in the study area using the buffering process, with a 50-metre radius around each house to discriminate the three main types of dwellings (isolated, scattered and clustered) and a 15-metre radius around each house belonging to clustered dwellings in order to discriminate the two sub classes of clustered dwellings (dense and very dense clustered dwellings). Three indicators were calculated: house density, mean area to be cleared per house (The French Forest Orientation Law of July 11 2001 makes clearing brush obligatory within a 50-metre perimeter around each house -although the local mayor can extend this limit to 100 metres -located at a distance of less than 200 metres from forest or scrubland) and mean perimeter per house to be protected by firefighters in the case of forest fire. Then, a statistical analysis was performed on the whole of polygons obtained by studying the entire study site. The aim was to assess whether there was any statistical difference, not only between the three main dwelling types based on spatial criteria such as housing density, area to be cleared, and perimeter to be protected, but also inside the dwelling cluster variable between dense and very dense clustered dwelling types. This calculation was made with Statgraphics software. Comparisons between the three populations corresponding to isolated, scattered and clustered dwellings were made as well as comparisons between the two sub populations of clustered dwellings, dense and very dense clustered dwellings.

Relationships between the distribution of fire ignition points and dwelling types

Regarding the distribution of fire ignition points in wildland-urban interfaces (WUIs) as a function of different types of dwelling, an indicator of fire ignition density was calculated using the fire ignition point database. This was defined as the ratio of the number of ignition points located within a given area to the surface of this area, and is expressed in the number of fires per 1,000 hectares. We considered different types of territory within the study area: isolated WUIs, scattered WUIs and clustered WUIs, and the remainder of the study area. The surface area of the different types was calculated and the number of fire ignition points located within each part was determined.

Results

Maps of dwelling types based on distance plus counting criteria using GIS

The process of creating buffer zones, based on both the spatial juxtaposition of houses and on the number of houses, enabled us to create a map of dwelling types using Arcgis®9.1 software. This is an efficient method for the characterization of dwelling types: isolated, scattered and clustered dwellings, and has the advantage of being easy to use. The method is based on distance and house-counting criteria using GIS and can thus be applied to large areas (fig. 2). Each house belongs to a single dwelling type, as we can see if we zoom in on the map (fig. 3). The repartition of polygons according to dwelling types is presented in table 1. Two sub classes (dense and very dense) were then defined in the clustered dwelling class, by buffering the houses included in the 72 polygons of this class. Fig. 4 compares the result of mapping dwellings with a photo corresponding to the extract of the mapped area. Maps of dwellings were produced easily and can be used in the first step in characterizing and mapping wildland-urban interfaces.

Dwelling types statistically different according to spatial criteria related to the risk of forest fire

Over the whole of the study area, the three populations corresponding to isolated dwellings (ID), scattered dwellings (SD) and clustered dwellings (CD) have been compared for each spatial criterion: mean housing density, mean area to be cleared and mean perimeter to be protected. The number of observations (corresponding to polygons) is 979 for the isolated dwelling population, 524 for the scattered dwelling population and 72 for the clustered dwelling population. A box plot comparing the three dwelling types for mean housing density (fig. 5.1), mean area to be cleared (fig. 5.3) and mean perimeter to be protected (fig. 5.5) shows that the type of dwelling has a significant effect on each response variable. Isolated dwellings have the highest mean area to be cleared and mean perimeter to be protected and the lowest mean housing density. In contrast, clustered dwellings have the lowest mean area to be cleared and mean perimeter to be protected and the highest mean housing density. Figures 5.1 to 5.6 present box plots, from left to right and from top to bottom, calculated for spatial criteria and different dwelling types (ID, SD CD and DCD, VDCD).

The two populations corresponding to the sub classes of clustered dwellings, dense clustered dwellings (DCD) and very dense clustered dwellings (VDCD) have been compared for each spatial criterion: mean housing density, mean area to be cleared and mean perimeter to be protected. The number of observations (polygons) is 411 for the dense clustered dwelling population and 142 for the very dense clustered dwelling population. As before, the box plots, comparing the two clustered dwelling types (dense and very dense) for mean housing density (fig. 5.2), for mean area to be cleared (fig. 5.4) or for mean perimeter to be protected (fig.5.6) show that the type of dwelling has a significant effect on response variable. Dense clustered dwellings have the highest mean area to be cleared and mean perimeter to be protected and the lowest mean housing density. In contrast, very dense clustered dwellings have the lowest mean area to be cleared and mean perimeter to be protected and the highest mean housing density. Fig. 5 (1)-( 6) Box plots (left to right, top to botton) Results confirmed that there were highly significant differences between the different types of dwellings. The non-parametric Kruskall-Wallis test was used to test both the difference between the three main classes of dwelling (isolated, scattered and clustered dwellings) and between the two sub classes of clustered dwelling (dense and very dense clustered dwelling). In each case the probability value is less than 0.05. As a result, there is a statistically significant difference for each criterion not only between the three main types of dwelling, but also between the two sub classes of clustered dwelling. Corresponding to highly statistically different types of dwelling, each of them has been characterized by their statistical values. Table 2 summarizes the main statistical values for housing density, area to be cleared and the perimeter to be protected for each type of dwelling. The actual values are presented in table 2 as: median-mean (standard deviation), minimum-maximum.

Table 2 Characteristics of dwelling types

Localisation of fire ignition points as a function of type of dwelling.

The mean fire ignition density value was calculated for the whole Aix-en-Provence study area for the period 1996-2007. This value reached 3.5 fire ignition points per 1,000 hectares. The same value was calculated specifically for wildland-urban interfaces, and proved to be 1.5 times higher in wildland-urban interfaces than elsewhere with 5.3 fire ignition points per 1,000 hectares (fig. 6). Our results show that fire frequency is very much higher in WUI areas than in natural areas. In addition, the fire ignition density was calculated for each of the dwelling types. Within dwelling types, we observed that the density value was two times higher for isolated dwellings than for clustered dwellings. Specifically, fire ignition density values increased greatly from clustered dwellings (4.2 fire ignition points per 1,000 ha), to scattered dwellings (5.2 fire ignition points per 1,000 ha) and finally to isolated dwellings (9.5 fire ignition points per 1,000 ha). This suggests that the spatial pattern of dwellings has a real impact on fire occurrence. Humans, and their spatial distribution, explain a part of the variability in the number of ignition points.

Discussion-Conclusions

The method developed to characterize dwelling types seems to be efficient, using GIS treatment based on the two following criteria: distance between houses and counting houses located close together. Moreover, this method is easy to use, can be transferred to land agencies or managers and can be applied to large areas on the scale of forest stands or several communities. The maps that can be produced using this method fulfil the needs of the main endusers concerned with forest fire prevention and forest fire-fighting. Indeed end-users such as forest and land planning managers or forest fire-fighters are interested in the location of wildland-urban interfaces that are directly related to different types of dwellings in contact with vegetation. Maps of different dwelling types can be interpreted for use in developing firefighting strategies or fire prevention measures.

Main criteria were defined: housing density, area to be cleared and the perimeter to be protected for each housing category. These criteria allow managers and fire-fighters to measure the effects of the dwelling structure on prevention measures or fire-fighting actions. It will be more difficult to protect isolated dwellings whose protection perimeter is almost ten times larger than the perimeters for clustered dwellings. It will be more difficult for a home owner to keep a large area around an isolated dwelling clear where this is three times more important than in a clustered dwelling. The mean density values that discriminate the three dwelling types proposed are approximately 100 houses per square kilometre for isolated dwellings, 160 houses per square kilometre for scattered dwellings, and 345 houses per square kilometre for clustered dwellings. These density values can be compared to thresholds defined in previous studies concerned with wildland-urban interface density values. [START_REF] Theobald | Land-use dynamics beyond the American urban fringe[END_REF] and [START_REF] Stewart | Characteristics and location of the Wildland-Urban-Interface in the United States[END_REF][START_REF] Stewart | Defining the Wildland-Urban Interface[END_REF] fixed two thresholds to differentiate wildland urban areas with a minimum housing density of 6 houses per square kilometre for WUI from urban areas with a housing density over 148 houses per square kilometre. Other thresholds were used by [START_REF] Hammer | Characterizing dynamic spatial and temporal residential density patterns from 1940-1990 across the North Central United States[END_REF], for whom the lowest-density rural clusters remain below 2 housing units per square kilometre, while the highest density category was more than 128 units per square kilometre. These thresholds can be combined with thresholds obtained with the density function [START_REF] Lampin | Dwellings characterization and mapping for forest fire risk prevention[END_REF]: less than 48 houses per square kilometre for isolated dwellings, around 48 to 160 houses per square kilometre for scattered dwellings, and more than 160 houses per square kilometre for clustered dwellings, where the highest values are relatively similar and the lowest values are different. This result is congruent with the low values obtained with the mathematical morphology approach [START_REF] Lampin | Dwellings characterization and mapping for forest fire risk prevention[END_REF]. Values obtained with our presented method are higher than the values specified in the methods. We can suppose both that the delimitation of the area that is the basis of the density calculation plays an important role in explaining this gap and that land cover and landscape are probably highly different.

The analysis of the spatial distribution of the fire ignition points shows a clear and positive relationship with the wildland-urban interface, as in Badia-Perpinya and [START_REF] Badia-Perpinya | Spatial distribution of ignitions in Mediterranean periurban and rural areas: the case of Catalonia[END_REF]. The results also show significant differences between the dwelling types. The spatial housing pattern influences the frequency of fire ignitions as in Cardille et al. (2001) and Pew and Larsen (2001). We point out that a lower housing density is connected with a higher number of fire ignition points. This is in accordance with the results of Syphard et al. (2007a), who used fine-resolution maps of WUI from housing density and land cover data. They showed that lowdensity housing explained significant variation in the number of fires.

In terms of fire management, these results should prove useful in identifying different dwelling types, which are an extremely important component of fire risk, in terms of the following two elements: hazard and vulnerability [START_REF] Jappiot | Ch.apitre 6 "Systèmes d'information géographique et modélisation dans le domaine de la prévention des incendies de forêt[END_REF][START_REF] Blanchi | Forest fire risk assessment and cartography. A methodological approach[END_REF]. Fire risk variables can be associated to two main groups: physical variables and human variables [START_REF] Amatulli | Assessing long-term fire ris kat local scale by means of decision tree technique[END_REF]. Dwelling density has only a small influence on fire risk when considered in the context of the number of possible socio-economic factors linked with the occurrence of wildfires [START_REF] Romero-Calcerrada | Gis analysis of spatial patterns of human-caused wildfire ignition risk in the SW of Madrid (central Spain)[END_REF]. According to [START_REF] Syphard | Human influence on California fire regimes[END_REF], ignition sources are linked to the spatial arrangement of dwellings and continuous vegetation. Other biophysical variables such as elevation and fuel types should be combined with the dwelling types in order to assess overall fire risk. Risk assessment could also be improved by including other socio-economic factors such as distance to the nearest road, population density, and others. In terms of public policy, introducing the risk of fire, and particularly the vulnerability of the territory using maps of dwelling types, is a way to make the inhabitants more aware of fire risk in WUI. This will globally decrease the risk of fire (i) by reducing fire propagation via biomass removal, and (ii) by reducing fire ignition probability by encouraging less carelessness. Accomplishing this goal is strictly related to the designation of suitable prevention messages and preventive measures which can be different according to WUI types. An important land use management implication of these results is to use this approach to measure urban sprawling, and to provide some suggestions to recommend where to encourage compact urban development. This is one of the current applications in France, namely to use this kind of map in Prevention Plans at the municipality level. Moreover, as [START_REF] Romero-Calcerrada | Gis analysis of spatial patterns of human-caused wildfire ignition risk in the SW of Madrid (central Spain)[END_REF] pointed out, these results could possibly be useful to fire managers, in order to concentrate their forces near the riskiest areas, and to thereby be able to attack a fire as quickly as possible. Finally, the qualification of dwelling types is one of the main elements of the characterization of wildland-urban interfaces, where this is crucial for the improvement of forest fire prevention and fire-fighting.

Fig. 1

 1 Fig.1 Location of the study area in the south-eastern France

Fig. 2

 2 Fig.2 Dwelling type map in the study area

Fig. 4

 4 Fig.4 Extract of mapped area corresponding to the photo

Fig. 6

 6 Fig.6 Fire ignition density value

Table 1

 1 Number of identified polygons in the study area

Acknowledgments

The work was partially financed by the French "Ministère de l'Ecologie, du Développement et de l'Aménagement Durables", "Conseil Régional Provence-Alpes-Côte-d'Azur". The authors wish to thank them. They also acknowledge J.C.Thill for reviewing the final manuscript and the anonymous peer reviewers for their comments, the firefighters of the "Service d'Incendie et de Secours du Var", and the foresters of the "Direction Territoriale Méditerranée de l'ONF" for their contribution.


Benoît Camenen

Cemagref Lyon, 3 bis quai Chauveau, CP 220, 69336 LYON cedex 09, FRANCE Email: benoit.camenen@cemagref.fr

The author discussed an important issue in mud floc knowledge, which is the uncertainties in the measurements of excess density and settling velocity. He pointed out the lack of accuracy in measuring instruments, and the statistical nature of the variables. However, the author did not discuss the limitations of the theories and semi-empirical formula he used. Some assumptions in the calculations need also to be discussed in more details; and some calculations are not mathematically correct.

The link between excess density and settling velocity

The fall velocity equation suggested by Winterwerp (1998, cf. Eq.5 in the original paper) relies on the Stokes law, i.e. it corresponds to the balance between the gravitational and drag force for a spherical non cohesive particle, and on an additional empirical equation, which includes the effect of the Reynolds number on the settling velocity (coefficient f (Re) = (1 + 0.15Re 0.687 ) -1 , see [START_REF] Graf | Hydraulics of Sediment Transport[END_REF]). The fall velocity may be written as follows :

w s = c f 18 ∆ρ g η D f 2 f (Re) (1)
where c f is a coefficient depending on the shape, roughness and porosity of the particles (c f ≈ α/β following Winterwerp developments). According to Eq.1, the standard deviation for the settling velocity is directly related to the standard deviation for the excess density ∆ρ = ρ f -ρ w . And there is no real need to introduce the fractal theory to explain the uncertainty of settling velocity. It only yields a new uncertainty as it introduces a new unknown parameter nf ... On the other hand, there are several additional terms that are missing in Eq.9 (original paper) such as the errors due to the coefficient c f , to the molecular viscosity of water, which is function of temperature and SPM concentration, and to the Reynolds number Re = w s D f /ν (ν is the kinematic viscosity of water), or more exactly to the function f (Re). The standard deviation σ ws should be written as follows :

σ ws = ∂w s ∂∆ρ 2 σ 2 ∆ρ + ∂w s ∂D f 2 σ 2 D f + ∂w s ∂c f 2 σ 2 c f + ∂w s ∂η 2 σ 2 η + ∂w s ∂f (Re) 2 σ 2 f (Re) (2)
If we follow the author assumptions, the three last terms are negligible, and so :

σ ws w s ≈ σ ∆ρ ∆ρ 2 + 4 σ D f D f 2 (3)
The author finally found as an average the relative σ D f /D f ≈ 30%, σ ∆ρ /∆ρ ≈ 10%, and then σ ws /w s ≈ 100%. If these results are coherent with Eq.3 (although overestimated, from Eq.3 σ ws /w s ≈ 60%), it clearly appears that the uncertaities on D f has a much larger influence on the uncertainties on w s than the uncertainties on ∆ρ. It should also be noted that Eq.1 (or Eq.5 in the original paper) is actually an implicit formula as Re is a function of w s ; thus, Eq.2 is not strictly correct (so are Eqs.9 and A4 to A8 in the original paper). Following the same idea, from a mathematical point of view, Eq.6 should be written :

σ ∆ρ = ∂∆ρ ∂M p 2 σ 2 Mp + ∂∆ρ ∂V f 2 σ 2 V f + ∂∆ρ ∂ρ p 2 σ 2 ρp + ∂∆ρ ∂ρ w 2 σ 2 ρw ( 4
)
as there is an uncertainty on ρ w (even if this latter is negligible compared to the other parameters).

On the estimation of the uncertainty

When one wants to discuss about uncertainties, it seems obvious that all uncertainties should be discussed. One cannot exclude systematic errors due to the lack of accuracy of the measuring instruments. Some attempt to estimate this error is needed; otherwise, the results are biased. For example, the volume concentration V f may typically be estimated as :

V f = c v πD f 3 6 (5)
where c v is a shape factor taking into account the fact that a floc is not spherical. Then, the standard deviation σ V f reduces to :

σ V f = ∂V f ∂D f 2 σ 2 D f + ∂V f ∂c v 2 σ 2 cv (6)
As the uncertainty on the floc size D f is much larger than the uncertainty on the floc sphericity and ∂V f /∂D f >> ∂V f /∂c v , σ cv may be neglected in Eq. 6.

For the settling velocity, an important uncertainty that needs to be discussed is σ c f . Based on a study on non-cohesive particle of various shape, [START_REF] Camenen | A simple and general formula for the settling velocity of particules in suspension[END_REF] showed that w s could be underestimated by 20% approximately. On the other hand, [START_REF] Johnson | Settling velocity of fractal aggregates[END_REF] showed that the settling velocity of fractal floc may be underestimated by a factor 5 because of the porosity of the floc. In the same way, the uncertainty on w s due to the empirical function of the Reynolds number σ f (Re) should be estimated :

σ f (Re) = ∂f (Re) ∂R e 2 σ 2 Re + σ 2 f (7)
where σ f is the standard deviation on the empirical equation f . In term of the presentation of the error, the author should clarify some points where results are getting unphysical. Defining the error on the primary particle size such as D p = 1.5 ± 3.5µm (cf. Tab.2) seems awkward (D p > 0). It may be better in that case to estimate an error on log D p . Similarly, it is incorrect to write nf = 3.23 ± 0.06 knowing that nf ≤ 3 by definition. However, it indicates that there is an error in the estimation of the uncertainties. Based on Eq.1 of the original paper (where ∝ may be replaced by = following [START_REF] Winterwerp | A simple model for turbulence induced flocculation of cohesive sediment[END_REF], a better estimation of nf would be :

nf = 3 + ln ρ f -ρw ρp-ρw ln D f Dp (8)
Using this equation, we ensure nf < 3 is satisfied as ρ f < ρ p and D f > D p . Eq.2 (original paper) yields too approximate results as D p cannot be assumed as an independent variable. From Eq. 8, the standard deviation σ nf is then :

σ nf = ∂nf ∂D f 2 σ 2 D f + ∂nf ∂D p 2 σ 2 Dp + ∂nf ∂ρ f 2 σ 2 ρ f + ∂nf ∂ρ p 2 σ 2 ρp + ∂nf ∂ρ w 2 σ 2 ρw (9)
Moreover, one should not confuse an uncertainty on the measurements of a variable and the physical distribution of this variable. It is somehow an error to discuss uncertainties on the diameter of flocs when measuring the size of a population of flocs.

Finally, as both settling velocity and floc density are the main unknown parameters for the cohesive sediments, it would be very interesting to combine the LISST 100C apparatus with some optical apparatus to determine the settling velocity as the INSSEV apparatus developed at Plymouth University [START_REF] Manning | The use of optics for the in situ determination of flocculated mud characteristics[END_REF]. Such an experiment would yield to significant improvements on the determination of the uncertainties for each method.

  Discussion of "Uncertainty of excess density and settling velocity of mud derived from in situ measurements" by M. Fettweis Estuarine Coastal & Shelf Science, Vol.78(2), June 2008 DOI: 10.1016/j.ecss.2008.01.007


François De Guio
email: francois.de-guio@cemagref.fr

Maja Musse

Hugues Benoit-Cattin

Tiphaine Lucas

Armel Davenel

Magnetic resonance imaging method based on magnetic susceptibility effects to estimate bubble size in alveolar products: Application to bread dough during proving

Keywords: magnetic susceptibility, MRI simulation, field inhomogeneity, gas bubble, gas cell, bread, dough, proving, fermentation

Magnetic resonance imaging has proven its potential application in bread dough and gas cell monitoring studies, and dynamic processes such as dough proving and baking can be monitored. However, undesirable magnetic susceptibility effects often affect quantification studies, especially at high fields. A new low field method is presented based on local assessment of porosity in spin-echo imaging, local characterization of signal loss in gradientecho imaging and prediction of relaxation times by simulation to estimate bubble radii in bread dough during proving. Maps of radii showed different regions of dough constituting networks which evolved during proving. Mean radius and bubble distribution were assessed during proving.

Introduction

Magnetic susceptibility effects are of increasing interest in Magnetic Resonance

Imaging (MRI). Traditionally, as field inhomogeneities induce geometry and intensity distortions in spin-echo (SE) imaging and signal loss in gradient-echo (GE) imaging [START_REF] Ludeke | Susceptibility artefacts in NMR imaging[END_REF][START_REF] Bakker | Simulation of susceptibility artifacts in 2D and 3D Fourier transform spin-echo and gradient-echo magnetic resonance imaging[END_REF], susceptibility artefacts have been considered undesirable and therefore subjected to correction [START_REF] Fernandez-Seara | Postprocessing technique to correct for background gradients in image-based R2* measurements[END_REF]. However, there are number of emerging techniques and applications using magnetic susceptibility effects and T 2 *-weighted imaging. The difference in magnetic susceptibility between deoxygenated and oxygenated blood, termed the BOLD effect, led to functional MRI which is widely used in cognitive neuroscience [START_REF] Ogawa | Functional brain mapping by blood oxygenation level-dependent contrast magnetic resonance imaging. A comparison of signal characteristics with a biophysical model[END_REF]. Susceptibility Weighted Imaging (SWI) uses the original phase image to enhance contrast between tissues with different susceptibilities [START_REF] Haacke | Susceptibility weighted imaging (SWI)[END_REF]. This technique has been applied to studies of brain tumors, trauma, vascular malformations and for quantification of brain iron [START_REF] Sehgal | Clinical applications of neuroimaging with susceptibilityweighted imaging[END_REF]. Superparamagnetic iron oxide (SPIO) particles create intense magnetic field distortions within and around cells, leading to irreversible signal dephasing in GE sequences. SPIO particles have been exploited as contrast agents for non-invasive cell tracking to determine their biodistribution in different organs [START_REF] Cunningham | Positive contrast magnetic resonance imaging of cells labeled with magnetic nanoparticles[END_REF][START_REF] Bos | In vivo MR imaging of intravascularly injected magnetically labeled mesenchymal stem cells in rat kidney and liver[END_REF][START_REF] Ittrich | In vivo magnetic resonance imaging of iron oxide-labeled, arterially-injected mesenchymal stem cells in kidneys of rats with acute ischemic kidney injury: Detection and monitoring at 3T[END_REF][START_REF] Mowat | Evaluating SPIO-labelled cell MR efficiency by three-dimensional quantitative T-2* MRI[END_REF]. Susceptibility studies have been applied to material characterization [START_REF] Beuf | Magnetic resonance imaging for the determination of magnetic susceptibility of materials[END_REF][START_REF] Peeters | Magnetic resonance imaging of microstructure transition in stainless steel[END_REF] and particle identification in industrial systems [START_REF] Robson | Identifying particles in industrial systems using MRI susceptibility artefacts[END_REF].

The potential of gas microbubbles as an in vivo intravascular susceptibility contrast agent for MRI has been demonstrated [START_REF] Wong | In vivo study of microbubbles as an MR susceptibility contrast agent[END_REF]. Moreover, differences in magnetic susceptibility (∆χ) between gas bubbles and a water environment create field inhomogeneities which induce intravoxel dephasing and associated signal loss characterized by a T 2 * shortening in GE imaging. However, the presence of air bubbles in several food products has not to date been the focus of studies of susceptibility effects. Understanding bubble mechanics is very important with a wide range of foamed food products. Bubble studies make possible the optimization of process design, contribute to the development of strategies for deaeration and are important in the area of texture and sensory analysis [START_REF] Niranjan | An Introduction to Bubble Mechanics in Foods[END_REF]. It can be anticipated that air or gas bubbles would be natural markers or contrast agents. Moreover, bubbles are part of a product and thus could be considered as a signature of the medium.

Proving flour dough is an essential stage in breadmaking. As the final character of most bakery products depends on the creation and control of gas bubble structures in the unbaked matrix, improving the understanding of the nucleation and growth of bubbles is of major interest [START_REF] Bellido | The bubble size distribution in wheat flour dough[END_REF]. MRI is a particularly suitable tool to study bread dough as it is non-invasive (essential for dough monitoring during proving as the alveolar structure is fragile and may collapse if intrusive measurement techniques are used), and provides a good trade-off between spatial and time resolution [START_REF] Bonny | Dynamic magnetic resonance microscopy of flour dough fermentation[END_REF][START_REF] Grenier | Assessment by MRI of local porosity in dough during proving. Theoretical considerations and experimental validation using a spinecho sequence[END_REF][START_REF] Ishida | Architecture of baked breads depicted by a magnetic resonance imaging[END_REF][START_REF] Van Duynhoven | Quantitative assessment of gas cell development during the proofing of dough by magnetic resonance imaging and image analysis[END_REF]. Dough samples can thus be quite large and imaged with high resolution during proving [START_REF] Van Duynhoven | Quantitative assessment of gas cell development during the proofing of dough by magnetic resonance imaging and image analysis[END_REF]. Even the entire baking process can be monitored with MRI [START_REF] Wagner | MRI study of bread baking: Experimental device and MRI signal analysis[END_REF]. Bonny et al. [START_REF] Bonny | Dynamic magnetic resonance microscopy of flour dough fermentation[END_REF] used high-field (9.4 T) magnetic resonance microscopy along with routine mathematical morphology to characterize the proving process noninvasively.

They faced critical susceptibility effects responsible for severe geometrical distortions due to the heterogeneous structure characterized by many interfaces. Different stages of fermentation could be identified, with a general increase in bubbles, but no quantitative assessment of bubble size could be made. Van Duynhoven et al. [START_REF] Van Duynhoven | Quantitative assessment of gas cell development during the proofing of dough by magnetic resonance imaging and image analysis[END_REF] illustrated the ability of MRI and image analysis to assess gas cell development in the growth of the dough during proving.

Basically, the size of each cell was determined by counting the number of pixels, thus presupposing high resolution and hence small fields of view. Moreover, this method does not account for susceptibility effects that are especially disrupting at 4.7 T and are responsible for geometry and intensity distortions in SE imaging. X-ray tomography has also proved its ability in analyzing bubble growth during bread making [START_REF] Babin | Fast X-ray tomography analysis of bubble growth and foam setting during breadmaking[END_REF]. However, this technique is only sensitive to density, and imaged samples have to be small (< 10 mm).

The aim of the study reported here was to characterize the alveolar structure of bread dough and its evolution during proving using susceptibility effects. Grenier et al. [START_REF] Grenier | Assessment by MRI of local porosity in dough during proving. Theoretical considerations and experimental validation using a spinecho sequence[END_REF] assessed local dough porosity (volume of gas per volume of dough) during proving. However, nothing can be inferred about bubble size for a given value of porosity, as different cell distributions can be obtained [START_REF] Bikard | Simulation of bread making process using a direct 3D numerical method at microscale: Analysis of foaming phase during proofing[END_REF][START_REF] Babin | X-ray tomography study of the cellular structure of extruded starches and its relations with expansion phenomenon and foam mechanical properties[END_REF]. We therefore developed an original method to estimate bubble size based on the local porosity assessed in SE imaging and on the susceptibility effects in GE imaging. MRI simulation was used to predict signal loss for virtual networks of gas bubbles embedded in a water matrix. A low magnetic field (BB 0 = 0.2 T) was used in order to avoid the severe geometric distortions encountered at high field in SE imaging and critical signal losses in GE imaging as the heterogeneous structures were characterized by many susceptibility interfaces [START_REF] Bonny | Dynamic magnetic resonance microscopy of flour dough fermentation[END_REF]. The dynamic side of bread proving requires acquisition to be quite fast, leading to a limited number of acquisitions, limited resolution and an appropriate hypothesis to compute maps of interest. Section 2 reports how relaxation times affected by susceptibility were predicted with simulation. Section 3 describes the materials and methods used with bread dough preparations, the MRI protocol, computation of the maps and the principle of the estimation algorithm. Section 4 presents the results on non-yeasted bread dough and yeasted dough during proving, followed by a discussion in section 5.

Prediction of relaxation times by simulation

The algorithm developed for estimation of bubble size was partly based on simulation results. This section introduces theoretical considerations and describes the simulation protocol. We demonstrate how relaxation times can be modeled at 0.2 T and what can happen at 1.5 T.

Theoretical considerations

We first verified that the static dephasing regime conditions for spherical magnetic perturbers of radius r as described by Yung et al. [START_REF] Yung | Empirical models of transverse relaxation for spherical magnetic perturbers[END_REF] were largely met, as 1/τ << δω with τ = r 2 /D and δω = γ∆χBB 0 /3, D being the translational diffusion coefficient of spins taken as the water diffusion coefficient and γ the proton gyromagnetic ratio. In fact, D is even smaller in bread dough.

Microscopic field inhomogeneities (i.e. magnetic field inhomogeneities over distances with orders of magnitude smaller than the voxel size) are responsible for irreversible signal decay (R 2 = 1/T 2 ) and mesoscopic field inhomogeneities (from perturbers smaller than the voxel size but greater than the diffusion length [START_REF] Fernandez-Seara | Postprocessing technique to correct for background gradients in image-based R2* measurements[END_REF][START_REF] Sukstanskii | Theory of FID NMR signal dephasing induced by mesoscopic magnetic field inhomogeneities in biological systems[END_REF]. Macroscopic field inhomogeneities are not present in T 2 * definition and induce non-exponential signal decay [START_REF] Fernandez-Seara | Postprocessing technique to correct for background gradients in image-based R2* measurements[END_REF][START_REF] De Guio | Signal decay due to susceptibility-induced intravoxel dephasing on multiple air-filled cylinders: MRI simulations and experiments[END_REF]. Macroscopic field inhomogeneities can typically arise from air inclusions or ferromagnetic objects. Finally, assuming a constant proton density across the voxel (ρ(r) = ρ), signal decay in a voxel of volume V is expressed in Eq. (1):

) contribute to R 2 # , the reversible portion of R 2 * (R 2 * = R 2 + R 2 # = 1/T 2 * )
) 2 ( sin ). 2 ( sin ). 2 ( sin . . ) ( * 2 z z y y x x t R voxel l tg c l tg c l tg c e V t S γ γ γ ρ - = (1)
where l x , l y , l z are the voxel dimension, and field inhomogeneities along one direction i are expressed with a linear gradient (ΔBB i = g i .i). In this study, we defined T 2 χ as the constant of the exponential curve fitting the initial part of the susceptibility-induced signal decay.

Microscopic, mesoscopic and macroscopic inhomogeneities were thus taken into account in the approximation of T 2 χ.

Simulation protocol

The purpose of the simulation was to investigate quantitatively the signal loss in GE MR images induced by networks of gas bubbles embedded in water. We used the SIMRI MRI simulator as described by Benoit-Cattin et al [START_REF] Benoit-Cattin | The SIMRI project: a versatile and interactive MRI simulator[END_REF]. In a previous study [START_REF] De Guio | Signal decay due to susceptibility-induced intravoxel dephasing on multiple air-filled cylinders: MRI simulations and experiments[END_REF], intravoxel modeling and associated signal decay were quantitatively assessed and experimentally validated in the case of a single well-resolved air-filled cylinder and in the case of a network of small interacting air-filled cylinders. An overview of the simulation framework is presented in Fig. 1. A 3D virtual object was defined by radius (r) and center to center distance between bubbles (ccd) expressed in pix obj (i.e. the pixel in the 2563 object space), difference in magnetic susceptibility (∆χ) and number of bubbles in each dimension (Nb). The additive property of the magnetic field, and the analytical formulation of field inhomogeneities created by one sphere [START_REF] Ludeke | Susceptibility artefacts in NMR imaging[END_REF], were combined to compute the relative field inhomogeneities (∆B/B 0 ).

From this map, SIMRI provides simulated MR images with susceptibility effects. At this stage, object parameters such as spin-lattice relaxation time T 1, spin-spin relaxation time T 2, and proton density ρ could be defined. NMR experiments (Minispec PC 120, Bruker SA, Wissembourg, France) on the bread dough showed several T 1 and T 2 components leading to multi-exponential decay, as previously described [START_REF] Ruan | Study of water in dough using nuclear magnetic resonance[END_REF][START_REF] Grenier | Assessment of NMR relaxation times and ice fraction in frozen dough[END_REF]. The first components were smaller the 256 object space), a 32x32 slice simulated image resulted in r = 0.5 pix pix ima (r = 0.125 pix ima , ccd = 0.375 pix ima with a 8x8 simulated image). For each object, GE images were simulated at different increasing TE from 4 to 12 ms. The mean gray level was then computed on regions of interest (ROI) enclosing several bubbles (visible on Fig. 1). This centred ROI made it possible to avoid boundary effects and to account for the interactions between several bubbles. Good representation of the medium was thus assured with such an ROI as bubbles were regularly spaced. T 2 χ was calculated as the constant of the exponential curve fitting the initial part of the mean gray level decay computed in the ROI. In fact, objects with large ccd gave quasi-exponential decay even for longer TE (12 ms). As ccd decreased, signal decay was hardly exponential because of high gradients g i in Eq. ( 1) and T 2 χ was assessed on the very first TE values (between 4 and 8 ms).

Evolution of T 2 χ as a function of bubble radius

Our aim was to link the object configuration with T 2 χ values so as to characterize the medium. The initial question was to find out whether small or large intravoxel bubbles engendered a similar signal loss for the same porosity (i.e. the gas volume fraction in the voxel). We therefore gathered all T 2 χ values relative to the different object configurations as a function of bubble radius at the different porosities shown in Fig. 2.a. The same curve shape was observed for each porosity. The first plateau indicated that a similar initial signal loss was measured for small radii, i.e. T 2 χ was only dependent on porosity. Then, a rapid increase in T 2 χ leveled off on a second plateau; note that the lower the porosity the smaller the radii corresponding to T 2 χ values as they started to increase. This was appropriate to our study as the initial mean porosity of bread dough was around 10% and overall porosity and bubbles increased during proving. One significant feature was that T 2 χ was sensitive to porosity and intravoxel structure according to bubble size.

Modeling of T 2 χ at 0.2 T

T 2 χ as a function of bubble radius was modeled with a sigmoid curve for each porosity value; a typical sigmoid function being defined in Eq. ( 2).

x e x f

λ - + = 1 1 ) ( (2)
As curves were not centered but shifted along the radius axis, we defined r break as the radius corresponding to the symmetrical point of the curve. As the amplitude between low and high radius values was not 1, K was introduced as multiplying factor. Finally, T 2 χ init was defined as T 2 χ offset, i.e. the value of T 2 χ for the smallest radius. Using these new variables and Eq. ( 2), T 2 χ as a function of radius is expressed in Eq. (3):

) ( 2 2 1 1 . ) ( break r r init e K T r T - - + + = λ χ χ (3)
For several porosity values, the parameters of Eq. (3) (T 2 χ init , K, r break , λ) were assessed to best fit the data according to the least mean square criterion (using Marquardt's algorithm in TableCurve2D® software, Jandel Scientific, version 5). Each parameter separately was then studied as a function of porosity. T 2 χ init and K (Fig. 3.a), r break and λ (Fig. 3.b) were drawn up and found to be closely related to porosity with second order polynomial functions (R 2 (T 2 χ init ) = 0.995, R 2 (K) = 0.988, R 2 (r break ) = 0.994 and R 2 (λ) = 0.993).

Finally, a complete description of T 2 χ as a function of radius and porosity (ε) is expressed in Eqs. [START_REF] Ogawa | Functional brain mapping by blood oxygenation level-dependent contrast magnetic resonance imaging. A comparison of signal characteristics with a biophysical model[END_REF][START_REF] Haacke | Susceptibility weighted imaging (SWI)[END_REF].

)) ( )( ( 2 2 1 1 ). ( ) ( ) , ( ε ε λ χ χ ε ε ε break r r init e K T r T - - + + = (4) With ( 5
+ + - = + - = - + - = + - = ε ε ε ε ε ε λ ε ε ε ε ε ε χ break init r K T
Note that the polynomial fitting on K as a function of ε was accurate and determined for porosity values greater than 10%. It was obvious that a homogeneous object without bubbles (ε = 0) would normally result in a zero multiplying factor (K). The same was valid for r break , whose limit for low porosities should be zero. Moreover, zero porosity in the T 2 χ init expression corresponded to the T 2 * value, i.e. the exponential signal decay constant representing the signal loss due to the mesoscopic field inhomogeneities.

Signal losses at 1.5 T

Although a low magnetic field was used, we felt that it could be of interest to predict signal losses in simulation at a higher field. The same simulations were therefore undertaken at 1.5 T and the results are summarized in Fig. 2.b. Susceptibility was obviously stronger as field inhomogeneities were proportional to the main magnetic field. When porosity was equal, T 2 χ init times were considerably shortened compared to the situation at B 0 = 0.2 T. Indeed, if ε = 0.12, simulations gave the followings: T 2 χ init = 12.65 ms at BB 0 = 0.2 T and T 2 χ init = 3.8 ms at B 0 B = 1.5 T. The value of BB 0 did not seem to impact on r break value. However, the K factor was significantly increased (1.35 at B 0 B = 0.2 T against 5.11 at BB 0 = 1.5 T if ε = 0.12). Signal decay was dramatically decreased for higher porosities (ε > 0.2), preventing GE studies at high fields.

We have shown in this section that T 2 χ relaxation time could be modeled with a sigmoid function dependent on porosity and radius of bubbles. We therefore used Eqs. [START_REF] Ogawa | Functional brain mapping by blood oxygenation level-dependent contrast magnetic resonance imaging. A comparison of signal characteristics with a biophysical model[END_REF][START_REF] Haacke | Susceptibility weighted imaging (SWI)[END_REF] in the next section to build an estimation algorithm to assess bubble size in bread dough.

Materials and methods

Sample preparation and experimental procedure

Bread dough samples were obtained by mixing 2000 g of wheat flour (Type 55, Moulins Soufflet Pantin), 1140 g of water, 40 g of salt, 20 g of improver and 100 g of yeast (optional) in a Moretti Forni Grain Spiry 8 dough mixer for 17 min at 100 rotations per minute. The advantages of using non-yeasted dough were to provide a stable object with regard to the acquisition time, with possible comparison of structure with yeasted dough at the initial time of proving. A fraction of gas (mainly air) is incorporated at the mixing stage and porosity at the end of mixing was estimated at 10 ± 2% depending on the mixer used for a

given recipe [START_REF] Cauvain | The Evolution of Bubble Structure in Bread Doughs and Its Effect on Bread Structure[END_REF], with a mean gas bubble size of approximately 0.05 mm to 0.3 mm [START_REF] Babin | Fast X-ray tomography analysis of bubble growth and foam setting during breadmaking[END_REF].

Dough temperature and water content were checked for the evaluation of the reproducibility between batches as both are known to affect the relaxation signal of dough [START_REF] Ruan | Study of water in dough using nuclear magnetic resonance[END_REF][START_REF] Grenier | Assessment of NMR relaxation times and ice fraction in frozen dough[END_REF]. The temperature was 24.5 ± 0.5°C, and water content was 45 ± 3 g of water per 100 g of dough. Cylindrical flasks (Ø = 50 mm, 70 mm height) were filled with 50 g of yeasted dough or 100 g of non-yeasted dough. A lid was placed over the flask to limit dehydration during measurement. Flasks were then placed in a tunnel within the magnet, equipped with thermal regulation set at 24.5°C. Internal tunnel temperature was monitored with thermocouples. MRI acquisitions were begun at approximately t= 7 min, t=0 referring to the end of mixing.

Cylindrical flasks were also filled with MnCl 2 solution, the concentration (1287.2 μM)

being adjusted to obtain a T 2 value close to that of bread dough, i.e. T 2 = 20 ms. SE images of the phantom were used to normalize dough MR images to correct for inherent magnet and coil spatial inhomogeneities.

MRI device and parameters

Porosity map

Assuming an initial uniform porosity at a reference time, the porosity of a voxel (i, j) can be estimated from its gray level during proving in SE imaging [START_REF] Grenier | Assessment by MRI of local porosity in dough during proving. Theoretical considerations and experimental validation using a spinecho sequence[END_REF]:

init init MGL j i SE j i ) , ( ). 1 ( 1 ) , ( ε ε - - = (6)
with ε init the initial overall porosity (ε init = 10%) at t init = 9.5 min, i.e. the time to obtain the first SE image, MGL init the corresponding mean gray level computed on a reference ROI, and SE(i,j) the gray level in the SE image at the location (i,j). We verified that the mean value computed in the porosity map was concordant with the overall porosity computed from the total dough volume measurement.

T 2 χ map

With GE 1 (i,j) and GE 2 (i,j) as the GE images at TE 1 and TE 2 , the T 2 χ map can be defined as follows:

) ) , ( ) , ( ln( ) ( ) , ( 2 1 1 2 2 j i GE j i GE TE TE j i T - = χ (7)
In fact, the exponential signal decay constant was estimated from two points only. Prior simulations showed that TE 1 = 4 ms and TE 2 = 12 ms were suitable for a 0.2 T magnetic field.

Bubble size estimation algorithm

The algorithm principle was to combine T 2 χ and porosity maps to estimate local bubble radii. Indeed, from Eq. ( 4) we can extract bubble radii (r) and compute the map of radii r(i,j) according to Eq. ( 8):

) 1 ) , ( ) , ( ) , ( ln( . ) , ( 1 ) , ( ) , ( 2 2 - - - = j i T j i T j i K j i j i r j i r init break χ χ λ (8)
r break , K, λ and T 2 χ init were determined for each pixel from the porosity map using Eq. ( 5).

T 2 χ was simulated from a set of several bubbles, i.e. signal loss was characterized in the center of the image to take into account the influence of the neighboring bubbles. This local approach based on several voxels should result in a coherent estimation algorithm, and therefore r(i,j) was not directly computed pixel-by-pixel but through a 3x3 averaging mask.

This allowed replacement of the pixel value by a local average around the pixel of interest.

Indeed, the non-linear property of the logarithm function made the subsequent use of an averaging on the map of radii impossible. Prior to computation of Eq. ( 8), we therefore applied the mask on r break , K, λ, T 2 χ and T 2 χ init maps.

Results

Non-yeasted dough

The porosity map of the non-yeasted dough was computed according to Eq. ( 6) and is represented in Fig. 4. The largest bubbles were easily detectable and the rest of the dough was quite uniform and dense. Mean porosity on a large reference ROI was found to be 12.6%. The corresponding mean value in the T 2 χ map computed on the same ROI was 12.9 ms. Compared with simulation results (T 2 χ = 13 ms and ε = 12.5 % in Fig. 2.a), these two values were in good agreement. The T 2 value of bread dough used for simulation and the general hypothesis of regularly spaced bubbles were justified to reach a good order of magnitude. The bubble size estimation algorithm was then applied to obtain the map of radii depicted in Fig. 5. The large bubbles were even more clearly revealed due to increased contrast compared to the neighboring environment. Fig. 5 also demonstrates localized regions of smaller bubbles which were hardly visible in the porosity map (Fig. 4). The map of radii depicted different regions corresponding to different bubble sizes. Smaller bubbles can be seen surrounded by larger bubbles, themselves surrounded by larger bubbles and so on. This radius evolution was spatially reproduced on the dough area, leading to different visible networks. The histogram representing the number of pixels according to their estimated radii extracted from the reference ROI is presented in Fig. 6. A Gaussian probability density function was found to fit well the experimental distribution of bubble size, as shown in Fig. 6 (σ = 0.062, μ = 0.37, R 2 = 0.991 with TableCurve2D® software, Jandel Scientific, version 5). The mean value computed in the reference ROI of the r(i,j) map was r = 0.38 pix ima . As detailed in the discussion section, estimated radius values were subjected to certain limits inherent to the method and values were not systematically significant. However, they can provide relative information about spatial differences (for a given protocol) and between different protocols of dough production.

Dough during proving

As in the previous experiment with the non-yeasted dough, r(i,j) was then computed for dough during proving. The three first maps at t 1 = 7 min, t 2 = 12.5 min and t 3 = 18 min after completing mixing are represented in Fig. 7 on the same scale as Fig. 5. Histograms were computed at t 1 , t 2 and t 3 , taking the same reference ROI occupying almost all the dough at t 1 (Fig. 8). As in the non-yeasted dough, distributions were Gaussian and shifted to the higher radius during proving, consistent with the growth of bubbles under desolubilization/vaporization of CO 2 , ethanol and water vapor. The average radius estimated from the algorithm as a function of the proving time is also presented in Fig. 9. Bubble size distribution was slightly shifted compared to non-yeasted dough, which may be attributed to the start of fermentation between the end of mixing and the first MRI acquisition (see Materials and Methods section). At t 2 , modification in the dough structure was clearly observed (Fig. 7.b). Groups of bubbles were growing and were identifiable from t 1 to t 3 and also for longer proving times. The overall networks, i.e. virtual limits between regions of the same radius, seemed to be almost unchanged in the horizontal direction during proving while bubbles underwent overall growth. From Fig. 7 and a complete set of maps of radii, dough evolution seemed to be related to the starting structure of the dough. Image processing methods on maps of radii would be useful for tracking purposes.

Discussion

In this section, we first discuss the simulation of T 2 χ and the behavior of the signal. We then explain the limitations of the method. Finally, results on non-yeasted dough and dough during proving will be analyzed in the light of this previous discussion.

Simulated evolution of T 2 χ

Simulation results concerning T 2 χ as a function of bubble radius (Fig. 2) were original and to our knowledge have never been reported in the literature. Several studies have been undertaken to study NMR signal dephasing due to the presence of mesoscopic field inhomogeneities in the static dephasing regime [START_REF] Yung | Empirical models of transverse relaxation for spherical magnetic perturbers[END_REF][START_REF] Sukstanskii | Theory of FID NMR signal dephasing induced by mesoscopic magnetic field inhomogeneities in biological systems[END_REF]. The relaxation constant was thus found to be independent of cylinder or bubble radii. From the free induction decay due to randomly distributed spherical particles, Yablonskiy [START_REF] Yablonskiy | Quantitation of intrinsic magnetic susceptibility-related effects in a tissue matrix. Phantom study[END_REF] derived an expression for the relaxation rate R 2 # defined in section 2.1:

0 # 2 . . . 3 9 4 R B χ γ ε π Δ = (9)
Eq. ( 9) explains T 2 χ shortening with increasing porosity (Fig. 3.a) or main magnetic field (Fig. 2.b), and the values are compatible with our present study. Indeed, if ε = 12.6%, Eq. ( 9) resulted in

T 2 * = 1 / (R 2 + R 2 # ) = 12.
2 ms and T 2 χ = 13 ms in our simulation.

We demonstrated three levels for T 2 χ observation. First, the plateau in Fig. 2.a indicated that, for a given porosity (or density), T 2 χ values are independent of bubble radii. This confirmed all the results encountered in mesoscopic scale studies. We showed that T 2 χ was dependent on radius in a certain range (macroscopic scale) and thus constitutes a sensitive indicator of the medium alveolar structure. The second plateau occurred for large bubbles compared to voxel size. T 2 χ was higher for larger bubbles because the signal from an air voxel could not decrease with TE due to the almost total absence of signal. In fact, heterogeneity was high in the ROI corresponding to the second plateau with complete air voxels, whereas for the first plateau each voxel contained small air bubbles. There were therefore more susceptibility interfaces for small bubbles, thus explaining the lower T 2 χ values. At a constant radius, an increase in porosity means a reduction in the center to center distance (ccd) and thus intravoxel dephasing is greater, as described in [START_REF] De Guio | Signal decay due to susceptibility-induced intravoxel dephasing on multiple air-filled cylinders: MRI simulations and experiments[END_REF] and in Eq. ( 1), explaining the decrease of T 2 χ with porosity.

Simulations at 1.5 T (Fig. 2.b) showed the potential of this method since there was a greater difference between T 2 χ init values for low porosities and an augmentation of the K factor. By using Eq. ( 8) to reveal bubbles, the estimated radius dynamic is thus enlarged and the algorithm would be more stable. Furthermore, a major gain in signal at high field would make better resolution possible. This is valuable for the study of the very small bubbles present in non-yeasted dough. Indeed, Bellido et al. reported a mean value of 100 μm in wheat flour dough [START_REF] Bellido | The bubble size distribution in wheat flour dough[END_REF].

Limitations of the method

Fig 2 .a illustrates some drawbacks of the method. T 2 χ is the same for small radii constituting the first plateau. Bubbles within this range cannot be discriminated. Fortunately, the process of dough proving with the growth of bubbles and the increase in local porosity was suited to this sensitive domain. Due to the shape of the T 2 χ curve (Fig. 2.a), computation of r(i,j) according to Eq. ( 8) worked well in the dynamic area between the two plateaux.

Although bread dough structure is highly heterogeneous (gaseous phase and dough films included), the model was quite accurate in terms of signal loss. A more accurate virtual description of the object will undoubtedly make possible the calculation of values of greater accuracy. In fact, the expansion of bubbles rapidly becomes heterogeneous and anisotropic (non-spherical and distorted bubbles) due to mechanisms of coalescence [START_REF] Bikard | Simulation of bread making process using a direct 3D numerical method at microscale: Analysis of foaming phase during proofing[END_REF]. However, using the estimation algorithm, it is possible to distinguish between bubbles according to their radius. At this point, r(i,j) values have to be considered as variables related to length rather than a physical variable. Babin et al. [START_REF] Babin | Fast X-ray tomography analysis of bubble growth and foam setting during breadmaking[END_REF] found dough structures were heterogeneous and dependent on the recipe, with bubble radii in the range 0.05-0.3 mm at the very beginning of dough proving. Maps of the radii in the present study were thus thought to be overestimated.

Non-yeasted dough and dough during proving

There is a real value in combining SE and GE sequences to assess local porosity and signal loss and thus deduce bubble radii. A porosity map (Fig. 4) provides important information (density is a criterion used by bakers to distinguish between bread recipes or to evaluate defects originating from the flour or the different stages in bread making) but a map of radii (Fig. 5) supplements the scientific understanding of evolution of the dough structure, from small to large bubbles. Additionally, the field of view available with MRI offers the possibility of assessing macroscopic heterogeneities (scale of a few millimeters to centimeters) in samples of realistic size comparable with industrial practices. Several image processing techniques based on this map could be developed to classify or segment bubbles.

The histogram presented in Fig. 6 shows a Gaussian distribution of bubbles. Bellido et al. [START_REF] Bellido | The bubble size distribution in wheat flour dough[END_REF] used microcomputed tomography on dough and found a log-normal distribution. These first results were encouraging for the study of non-yeasted bread dough.

Bikard et al. [START_REF] Bikard | Simulation of bread making process using a direct 3D numerical method at microscale: Analysis of foaming phase during proofing[END_REF] analyzed the foaming phase during proving using 3D simulation. They characterized the influence of parameters such as dough viscosity, kinetics or initial number of bubbles on the evolution of average bubble radius over time. Fig. 9 obtained from maps of radii at different times of proving is in good agreement with their results while looking at curve shapes. From a quantitative point of view, it was hard to compare as values were highly dependent on the above parameters that we could not quantify in our experiments.

Using digital image analysis, Rouille et al. [START_REF] Rouillé | Effect of flour minor components on bubble growth in bread dough during proofing assessed by magnetic resonance imaging[END_REF] showed that reduction in the number of small bubbles (Ø < 1 mm) was proportionally balanced by an increase in the number of larger bubbles ( 1 mm < Ø < 2 mm) during proving. This is also obvious from Figs. 78.

Image analysis-based methods [START_REF] Bonny | Dynamic magnetic resonance microscopy of flour dough fermentation[END_REF][START_REF] Van Duynhoven | Quantitative assessment of gas cell development during the proofing of dough by magnetic resonance imaging and image analysis[END_REF] showed limitations due to small sample size (magnetic resonance microscopy) and long imaging times. Even for small pixel size (115 μm), resolution is not sufficient to detect small bubbles at the beginning of the proving time [START_REF] Rouillé | Effect of flour minor components on bubble growth in bread dough during proofing assessed by magnetic resonance imaging[END_REF]. Susceptibility-based methods are able to reveal information at a lower scale due to the expansion of field inhomogeneities and associated signal loss.

Conclusion

Because they arise from object-dependent field inhomogeneities, magnetic susceptibility effects can be a source of quantitative data to characterize alveolar products such as bread dough. While classic image analysis techniques are limited by image or temporal resolution, this new method combines local porosity, local signal loss and simulation predictions to assess bubble radii which can be smaller or the same order of magnitude than the voxel.

However, fully quantitative accuracy of the method would suppose comparisons between different techniques or the use of a test object with known geometry and susceptibility.

Moreover, further simulation studies with more realistic object geometries to mimic bread dough would make possible more quantitative results. Growth and distribution of bubbles were observed during proving. Maps of radii gave information about dough structure and evolution. The principle of the estimation algorithm was shown at 0.2 T and simulations predicted greater distinction between bubble sizes at 1.5 T, especially for low porosities which are also encountered at key stages in breadmaking. Second order polynomial relationship with porosity was found for each simulation parameter.

  than the echo times (TE) used in this study. As a single value for relaxation times is necessary in simulation, weighted averaging of the relaxation times of the different components was performed, resulting in T 1 = 100 ms and T 2 = 20 ms. Lodi et al.[START_REF] Lodi | Characterization of water distribution in bread during storage using magnetic resonance imaging[END_REF] assessed T 2 values from spin-echo images by fitting time-series points pixel-by-pixel, resulting in T 2 maps of soy bread (mean values around 18 ms) consistent with this averaging. Simulation parameters such as main magnetic field (BB 0 = 0.2 T), pixel bandwidth (BW = 279 Hz/pixel) and repetition time (TR = 300 ms) were set to match real experimental conditions. By defining the number of pixels constituting the simulated image, we could define the relative proportion of the bubble radius in the simulated image pixel (pix ima ). For instance, considering the high-resolution object used as an example in Fig.1(r = 4 pix obj , ccd = 12 pix obj , ∆χ = -9.05 ppm, Nb = 17 in

A 0. 2 T

 2 electromagnet scanner in open configuration (Magnetom Open, Siemens, Erlangen, Germany) equipped with a head coil was used to image the bread dough. Double GE sequences (TE 1 = 4 ms, TE 2 = 12 ms) were performed to assess signal loss and T 2 χ during dough proving. SE sequences (TE = 8 ms) were used to obtain local porosity. For all sequences, slice thickness (ST), bandwidth (BW), field of view (FOV), matrix size (N), number of accumulations (Acc) and repetition time (TR) were set as follows: ST = 5 mm, BW = 279 Hz/pixel, FOV = 128 x 128 mm 2 , N = 128x128, Acc = 4, TR = 300 ms. Each sequence lasted 2 min 33 s, and GE and SE sequences were alternated.

Figure 1 :

 1 Figure 1: Overview of simulation framework. Field inhomogeneities map was dependent on

Figure 2 :B= 1 . 5 T

 215 Figure 2: T 2 χ value as a function of bubble radii at different porosities with sigmoid fitting

Figure 3 :

 3 Figure 3: a) K and T 2 χ init as functions of porosity. b) r break and λ as functions of porosity.

Figure 4 :

 4 Figure 4: Porosity map of non-yeasted bread dough. Large bubbles can be distinguished as

Figure 5 :

 5 Figure 5: Map of estimated radii of bubbles in non-yeasted bread dough. Contrast

Figure 6 :

 6 Figure 6: Histogram with Gaussian distribution (σ = 0.062, μ = 0.37, R 2 = 0.991) of the non-

Figure 7 :

 7 Figure 7: Map of estimated radii of bubbles in bread dough during proving. a) At t = 7 min b)

Figure 8 :

 8 Figure 8: Histograms extracted from maps of radii showing distribution of bubble sizes at

Figure 9 :

 9 Figure 9: Average estimated radius as a function of proving time.

Fig2a

  Fig2a

Fig.2b

  Fig.2b

Fig.3a

  Fig.3a

Fig.3b

  Fig.3b

Fig. 9

 9 Fig.9

ima , ccd = 1.5

Acknowledgments

The authors would like to thank D. Le Ray for dough preparations and the Regional Council of Brittany for financial support.