Anne-Sophie Bonnet-Lebrun
email: anne-sophie.bonnet-lebrun@normale.fr

Andrea Manica

Ana S L Rodrigues

Effects of urbanisation on bird migration

Keywords: anthropogenic effects, human population, partial migration, North American passerines, residency

Nearly one in five bird species is migratory, but not all individuals within a migratory species necessarily migrate: in partially migratory species, some do and some do not. Such within-species variability provides a natural experiment for investigating the mechanisms driving bird migration. Previous studies at the species level suggest that migrating provides a way to escape harsh winters, and to secure an increased access to resources, particularly important during the breeding season. Urbanization, by altering local temperatures ('heat island' effect) and resource availability (e.g. through garbage or garden feeders) can buffer the effects of winter harshness and modify breedingseason resource availability, potentially affecting individual migratory strategies. Here, we use ringing data from twelve North American partially migratory bird species to investigate the effects of natural environmental conditions (winter temperature, breeding season resource surplus) and urbanization on the propensity of individuals to migrate. We find strong support for the hypothesis that individuals migrate to avoid harsh winters, with, for eleven species, significantly higher probabilities of residency in areas with milder winters. We also found (significant for five species) that resource surplus in the breeding season reduces the propensity to migrate. Finally, urbanization increased the likelihood that individuals remain year-round in their ranges, avoiding to migrate away from their breeding range (four species) or their wintering areas (eight species), after controlling for climate and resources. Our results thus indicate that bird migratory strategies will respond to global changein climate and land useand indeed are already doing so.

Introduction

Nearly one in five of the world's 10,000 bird species migrate seasonally between breeding and nonbreeding ranges [START_REF] Kirby | Key conservation issues for migratory land-and waterbird species on the world's major flyways[END_REF], a global-scale ecological readjustment [START_REF] Moreau | The Place of Africa in the Palaearctic Migration System[END_REF]) that radically changes the composition and diversity of bird communities across wide areas of the planet [START_REF] Somveille | Mapping global diversity patterns for migratory birds[END_REF]. In practice, though, the migratory movements are not of species but of individuals. Furthermore, it is not a fixed species trait: in many species, some individuals migrate while others remain in the same region as year-round residents. This includes variation among populations whereby some populations are resident and some migratorybut also within populationsin which only a fraction of the individuals found in a given region migrate. The term 'partial migration' is often applied to the latter [START_REF] Chapman | The ecology and evolution of partial migration[END_REF], but here we use it more broadly to refer to within-species variation in migratory-versus-resident behaviour.

Partial migration is widespread among animal taxa (e.g. ungulates, [START_REF] Hebblewhite | Demographic balancing of migrant and resident elk in a partially migratory population through forage-predation tradeoffs[END_REF]fish, Chapman et al. 2012;insects, Dällenbach et al. 2018) and very common in bird species. For example, in a review of Australian land birds, [START_REF] Chan | Partial migration in Australian landbirds: a review[END_REF] found that among 155 non-passerine and 317 passerines species studied, respectively 44% of and 32% are partial migrants. Partially migratory species provide opportunities for testing hypotheses on the ecological and evolutionary processes underpinning migration itself, by allowing for multiple replicates (i.e., the individuals) while controlling for the wider variation in ecological traits (e.g. trophic level, body size, habitat preferences) observed across species (see [START_REF] Chapman | The ecology and evolution of partial migration[END_REF] for a review). Furthermore, given that individuals and populations can react faster to environmental changes than entire species, monitoring the responses of partially migratory species to anthropogenic activities can provide early insights into how species respond to global change [START_REF] Pulido | Current selection for lower migratory activity will drive the evolution of residency in a migratory bird population[END_REF][START_REF] Podhrázský | Central European Greylag Geese Anser anser show a shortening of migration distance and earlier spring arrival over 60 years[END_REF].

Conceptually, an individual bird in any given (breeding or non-breeding) location, can either remain in that same location as a resident, or move elsewhere for the following season. Recent studies at the species' level suggest that migration is largely driven by energetic trade-offs, with species migrating when the benefits derived from spending parts of the year in different areas exceed the costs of migration between them [START_REF] Hurlbert | The Effect of Energy and Seasonality on Avian Species Richness and Community Composition[END_REF][START_REF] Dalby | Seasonality drives global-scale diversity patterns in waterfowl (Anseriformes) via temporal niche exploitation[END_REF][START_REF] Somveille | Why do birds migrate? A macroecological perspective[END_REF]Somveille et al. , 2018a, b), b). More specifically, these studies suggest that the main driver of migration out of the breeding range is winter harshness, either because of the thermoregulation costs of low temperatures, or because of a reduction in resources during winter, or both [START_REF] Herrera | On the breeding distribution pattern of European migrant birds: MacArthur's theme reexamined[END_REF][START_REF] Lemoine | Potential Impact of Global Climate Change on Species Richness of Long-Distance Migrants[END_REF][START_REF] Carnicer | Geographic Differences between Functional Groups in Patterns of Bird Species Richness in North America[END_REF][START_REF] Schaefer | Impact of climate change on migratory birds: community reassembly versus adaptation[END_REF][START_REF] Somveille | Why do birds migrate? A macroecological perspective[END_REF]. At the individual level, this should translate into a higher probability of individuals to migrate out of their breeding locations if these locations face harsh winters, and conversely into a higher probability of remaining as resident in locations with mild winters. Accordingly, progressively milder winters linked to climate change have been proposed as a reason for increased rates at which some migratory species are being observed year-round in their breeding areas (e.g. Great Crested Grebes Podiceps cristatus in the Netherlands, [START_REF] Adriaensen | Ringing recoveries and the increase in numbers of European Great Crested Grebes Podiceps cristatus[END_REF]blackbirds Turdus merula in Europe, Berthold 1993;[START_REF] Main | Obligate and facultative partial migration in the Blackbird (Turdus merula) and the Greenfinch (Carduelis chloris): uses and limitations of ringing data[END_REF], likely reflecting an increase in the fraction of sedentary individuals in partially migratory populations (e.g. as observed for blackbirds in the Netherlands, [START_REF] Vliet | Changes in migration behaviour of Blackbirds Turdus merula from the Netherlands[END_REF]; and Denmark, [START_REF] Kristensen | Urbanization effects on migration behaviour of short-distance migratory birds[END_REF]. A relationship between winter harshness and migratory propensity can also explain observations that some previously non-migratory species whose ranges recently expanded into higher latitudes have become migratory there (e.g., European serins Serinus serinus spreading into northern Europe, [START_REF] Berthold | A comprehensive theory for the evolution, control and adaptability of avian migration[END_REF]; House Finch Carpodacus mexicanus in eastern North America, [START_REF] Able | Rapid 'evolution' of migratory behaviour in the introduced house finch of eastern North America[END_REF]).

The main drive to leave the wintering grounds, on the other hand, seems to be to obtain better access to resources during the breeding season. Indeed, previous studies at the species' level found that the resources available during the breeding seasonor, more precisely, the surplus in relation to the resources in the non-breeding seasonwas the main predictor of the number of breeding migratory birds found in any given region, a pattern explained in terms of reduced competition with resident birds [START_REF] Dalby | Seasonality drives global-scale diversity patterns in waterfowl (Anseriformes) via temporal niche exploitation[END_REF][START_REF] Somveille | Why do birds migrate? A macroecological perspective[END_REF]. At the individual level, one thus expects that the individuals wintering in areas with high levels of resource surplus in the breeding season have a higher probability of remaining as residents.

Given that anthropogenic activities are affecting both the climate (IPCC, 2014) and the distribution of natural resources [START_REF] Haberl | Human Appropriation of Net Primary Production: Patterns, Trends, and Planetary Boundaries[END_REF], it is not surprising that there is mounting evidence of human effects on bird migration [START_REF] Visser | Climate change leads to decreasing bird migration distances[END_REF][START_REF] Plummer | Is supplementary feeding in gardens a driver of evolutionary change in a migratory bird species?[END_REF][START_REF] Greig | Winter range expansion of a hummingbird is associated with urbanization and supplementary feeding[END_REF]. One such major effect is by altering land cover, with urbanization in particular creating habitats that are in many ways radically distinct from natural ones. The resulting changes to local environmental conditions may plausibly affect the migratory decisions of urban birds. First, urban areas are often 'heat islands', multiple degrees warmer than surrounding areas [START_REF] Collier | The impact of urban areas on weather[END_REF], which may buffer birds against harsh winter temperatures [START_REF] Shochat | From patterns to emerging processes in mechanistic urban ecology[END_REF]) and thus increase their probability of remaining over winter. Second, urban areas can provide particular resources, for example through bird feeders, garbage, and garden flowers and fruits (e.g. Robb et al. 2008;[START_REF] Greig | Winter range expansion of a hummingbird is associated with urbanization and supplementary feeding[END_REF]. These may either buffer populations against the low resources of winter, and thus increase the probability of breeding birds remaining over winter, or provide resources in the breeding season that increase the odds that wintering birds remain to breed. Overall, urbanization is thus expected to increase the propensity to residency among migratory birds.

Previous studies have already found evidence for such an effect, across a range of species. For example, evidence from ring recoveries [START_REF] Kristensen | Urbanization effects on migration behaviour of short-distance migratory birds[END_REF], stable isotopes [START_REF] Evans | Colonisation of urban environments is associated with reduced migratory behaviour, facilitating divergence from ancestral populations[END_REF], and physiological and behavioural studies [START_REF] Partecke | Increased Sedentariness in European Blackbirds Following Urbanization: A Consequence of Local Adaptation[END_REF] of European blackbirds (Turdus merula) indicates that urban birds have a lower tendency to migrate in relation to their nonurban counterparts. This effect is stronger at the northern part of the range [START_REF] Evans | Colonisation of urban environments is associated with reduced migratory behaviour, facilitating divergence from ancestral populations[END_REF], and seems to have played an important role in the relatively recent northwards expansion of the area in which the species is found year-round [START_REF] Møller | Loss of migration and urbanization in birds: a case study of the blackbird (Turdus merula)[END_REF]. Similarly, analysis of European robins (Erithacus rubecula) across different habitats near Antwerp, Belgium, found only 30% of colourmarked breeding birds in a woodland plot remained locally over winter, whereas most birds (including nearly all males) did so in urban garden and urban park plots [START_REF] Dhondt | Population Dynamics and Partial Migration of the European Robin (Erithacus rubecula) in Different Habitats[END_REF]. At a multispecies level, a study of wintering bird communities in Poland found higher abundances in urban compared to rural areas [START_REF] Tryjanowski | Winter Bird Assemblages in Rural and Urban Environments: A National Survey[END_REF]. Some studies investigated the mechanisms behind changes in migratory behaviour. Among these, some highlighted the role of increasing resources on bird migration. For example, bird feeders were found to increase the frequency of wintering Eurasian blackcaps (Sylvia atricapilla) in urban areas in Britain [START_REF] Plummer | Is supplementary feeding in gardens a driver of evolutionary change in a migratory bird species?[END_REF] and of Carolina wrens (Thryothorus ludovicianus) staying after winter at the northern edge of their range in Michigan [START_REF] Job | Wrens on the edge: feeders predict Carolina wren Thryothorus ludovicianus abundance at the northern edge of their range[END_REF]. Also, the use of nectar feeders seems to play a key role in the northwards expansion of wintering Anna's hummingbirds (Calypte anna) in the United States [START_REF] Greig | Winter range expansion of a hummingbird is associated with urbanization and supplementary feeding[END_REF], whereas overwintering birds in a population of mainly breeding migrant white storks (Ciconia ciconia) rely strongly on landfills [START_REF] Gilbert | Are white storks addicted to junk food? Impacts of landfill use on the movement and behaviour of resident white storks (Ciconia ciconia) from a partially migratory population[END_REF]. Some studies' results support the hypothesis that urban areas buffer migratory birds from winter harshness. For example, a study of ringing records of European blackbirds in the Netherlands found that the fraction of birds migrating away from their breeding areas in the winter was inversely correlated with temperature in rural areas, but found no noticeable effect in urban areas [START_REF] Vliet | Changes in migration behaviour of Blackbirds Turdus merula from the Netherlands[END_REF]). The urban 'heat island effect' was also considered a factor explaining the above-mentioned winter range expansion of Anna's hummingbirds in the United States, with the interaction between January minimum temperatures and housing densities being a significant predictor of winter presence [START_REF] Greig | Winter range expansion of a hummingbird is associated with urbanization and supplementary feeding[END_REF].

Most of previous studies have however analysed the effects of urbanization on bird migration by focusing on small parts of the range (e.g. [START_REF] Dhondt | Population Dynamics and Partial Migration of the European Robin (Erithacus rubecula) in Different Habitats[END_REF][START_REF] Gilbert | Are white storks addicted to junk food? Impacts of landfill use on the movement and behaviour of resident white storks (Ciconia ciconia) from a partially migratory population[END_REF][START_REF] Greig | Winter range expansion of a hummingbird is associated with urbanization and supplementary feeding[END_REF], and typically only at a specific season (usually winter). Additionally, measures of migratory propensity are oftentimes crude (e.g. contrasts between local bird densities, [START_REF] Tryjanowski | Winter Bird Assemblages in Rural and Urban Environments: A National Survey[END_REF]; fraction of colour-marked birds, [START_REF] Dhondt | Population Dynamics and Partial Migration of the European Robin (Erithacus rubecula) in Different Habitats[END_REF]; stable isotopes sensitive to only very large differences in migratory distance, [START_REF] Evans | Colonisation of urban environments is associated with reduced migratory behaviour, facilitating divergence from ancestral populations[END_REF], and can thus mask important differences in migratory behaviour (e.g. local post-breeding dispersion versus long-term migration).

Here, we take advantage of a continental-scale bird ringing scheme to test the hypothesis that urbanization is affecting bird migrations. Ringing data allow us to identify the precise seasonal locations of individual birds, and therefore unambiguously distinguish residency or short-distance dispersal from long-distance movements. We focus on ten North American partially migratory species for which we were able to obtain ringing recoveries data that cover their entire range. Furthermore, we investigate drivers in the migratory propensity of individuals in both seasons, breeding and nonbreeding, when accounting both for natural factors (winter harshness and resource availability) and for anthropogenic factors (urbanization). Specifically, we test four hypotheses: that propensity to remain year-round in the breeding grounds is higher if (1a) winters are milder and (1b) urbanization levels stronger; and that propensity to remain year-round in the wintering grounds is higher if (2a) resource surplus is higher and (2b) urbanization is higher.

Data and Methods

Individual bird data

Species' data come from the North American Bird Banding Programme (NABBP), run by the Bird Banding Laboratory of the USGS Patuxent Wildlife Research Center and the Bird Banding Office of the Canadian Wildlife Service (USGS Bird Banding Laboratory 2016). The programme compiles records of birds captured and marked with a uniquely numbered band or ring, as well as any subsequent recoveries. Ringing localities are recorded in 10-arc minute blocks, which, at 40°N for example, corresponds to 14.2km in longitude and 18.5km in latitude.

We focused on migratory species whose entire life cycle is well covered by the NABBP, i.e., whose breeding and non-breeding ranges (in the Western Hemisphere) fall mainly within the United States and southern Canada. We then searched within each species for individuals seen alive at least twice and in opposite seasons, i.e., one record in the breeding season (May to July), another in the nonbreeding season (December to February). Whenever an individual was recorded more than once in a season, we retained only the first record. Using the great circle distance between the breeding and nonbreeding locations, we classified each individual into either 'resident' (< 20km) or 'migrant' (>100km; following [START_REF] Fiedler | S19-3 Bird ringing as a means for studying the evolution of avian migrationpotentials, limits and prospects[END_REF]; [START_REF] Brown | Band recoveries reveal alternative migration strategies in American Robins[END_REF]). We focused on twelve species for which we were able to obtain a reasonable number of individuals (≥ 50) including a mix of resident and migrant individuals: American Goldfinch (Spinus tristis), American Robin (Turdus migratorius), Blue Jay (Cyanocitta cristata), Brown-headed Cowbird (Molothrus ater), Common Grackle (Quiscalus quiscula), European Starling (Sturnus vulgaris), Evening Grosbeak (Coccothraustes vespertinus), House Finch (Haemorhous mexicanus), Pine Siskin (Carduelis pinus), Purple Finch (Haemorhous purpureus), Red-winged Blackbird (Agelaius phoeniceus) and White-throated Sparrow (Zonotrichia albicolis). The eastern House Finch population has been recently introduced and is mostly resident [START_REF] Able | Rapid 'evolution' of migratory behaviour in the introduced house finch of eastern North America[END_REF], so we exclude it from the analysis, by removing all House Finch records east of 95°W. Total number of birds and number of residents are presented in Fig. 1). Records span the 1920-2016 period (Appendix A). Ringing/recovery effort, and thus NABBP records, are spatially clustered in regions of higher human density and those with more nature-friendly habitants, but we have no reason to expect ringing/recovery effort to affect the propensity of individuals to migrate.

Environmental data

Winter harshness [START_REF] Somveille | Why do birds migrate? A macroecological perspective[END_REF] found that the diversity of bird species that overwinter in their breeding locations was inversely related to both the winter temperature and to winter resources (measured by mean winter NDVI), with both variables being highly correlated and the former being the best predictor. Accordingly, we focused on winter temperature as an indicator of harshness of conditions during the non-breeding season. For any given location, we measured winter harshness as the mean of the monthly average temperature values from December to February within a 10km buffer around the location. To avoid focusing on an unrepresentative year, we used the Worldclim dataset (Worldclim database at resolution 30''; [START_REF] Hijmans | Very high resolution interpolated climate surfaces for global land areas[END_REF], averaging values across all available years .

Breeding resource surplus

Following [START_REF] Hurlbert | The Effect of Energy and Seasonality on Avian Species Richness and Community Composition[END_REF] and [START_REF] Somveille | Why do birds migrate? A macroecological perspective[END_REF]Somveille et al. ( , 2018a, b), b), we used values of mean monthly Normalized Difference Vegetation Index (NDVI) as a general indicator of resources (food, nesting sites and roosting sites). For any given location, we first calculated the mean monthly NDVI in the breeding season (May to July) and in the wintering season (December to February) within a buffer of 10km around the location, and then measured the resource surplus in the breeding season (ΔNDVI) as the difference between the former and the latter. Assuming (in a simple way) that resident species use the same resources year-round, this surplus is a measure of the resources available to migratory species (Somveille et al. 2018a, b). We obtained mean monthly NDVI values from NASA's Earth Observatory (2016; resolution 0.1˚), again averaging across all available years (May 2000 -February 2016) to create seasonal means.

Urbanization levels

We considered local human population density as a proxy for the level of urbanization. The median ringing/re-sighting year for the pooled species data was 1956 (interquartile: 1946-1968; Appendix A), so we used the closest spatially explicit data on reconstructed population density across North America that we were able to obtain (1970; Center for International Earth Science Information Network -CIESIN -Columbia University, 2017; resolution 30''). We extracted the mean population density in a buffer of 10km around each individual location, and log transformed it using log(x+1) (Fig. 1C,1E). A log transformation improved the distribution of this variable, which spans several orders of magnitude.

Testing of hypothesis

We tested four hypotheses: that the probability of an individual remaining year-round at its breeding location is (1) higher if winters are milder (i.e., higher local temperatures during the wintering season) and (2) higher if urbanization levels are stronger (i.e., higher local human density); and that the probability of an individual remaining year-round at its wintering location is (3) higher if breeding resource surplus is higher (i.e., higher local ΔNDVI) and (4) higher if urbanization levels are stronger (i.e., higher local human density).

We tested these hypotheses for each species by fitting a binomial Generalised Linear Model (logit link), with the resident (1) vs. migrant (0) status of each individual bird as a response and local environmental conditions (natural: winter temperature, ΔNDVI; anthropogenic: human density) as predictors. We tested hypotheses 1 and 2 by focusing solely on the breeding locations, using as predictors local winter temperature and local (log-transformed) human density, and hypotheses 3 and 4 by focusing on the wintering locations, using local breeding resource surplus (ΔNDVI) and local human density as predictors. Note that we did not model the probability of birds being present at a given location, but rather their probability of remaining as residents at a given location, knowing that they were present in a given season.

We standardised variables prior to modelling to allow for the comparison of estimated coefficients, using the following formula: (xmean(x))/sd(x). We used a backwards stepwise selection procedure using AIC to select the best model, and a sequential Bonferroni correction to deal with the large number of repeated tests across species. All analyses were done in R.3.4.2 (R Core Team 2017). Codes are provided in Appendix S4.

Results

Probability of remaining at the breeding location over winter

We found for eleven out of twelve species a significant positive relationship between the probability that individuals overwinter in their breeding locations and local winter temperature (Fig. 1B,2A, Appendix B). Human density was significantly positively related with the probability that individuals remain at their breeding locations overwinter for four out of twelve species: American Goldfinch, European Starling, Evening Grosbeak and Purple finch (Fig. 1C,2B). The effect was positive for five other species, and negative for three but not significantly so after sequential Bonferroni correction (Fig. 2B,Appendix B).

Probability of remaining at the winter location during the breeding season

For five out of twelve species -American Goldfinch, American Robin, Common Grackle, European Starling, and Evening Grosbeakwe found a significant positive effect of the surplus in resources during the breeding season on the probability that individuals remain at their wintering location into the breeding season (Fig. 1D,2C,Appendix C). This effect was positive but not significant for six other species, and negative but not significant for the remaining species (Red-winged Blackbird) (Fig. 2C,Appendix C). We found a significant effect of local human density on the probability that individuals remain at their wintering location during the breeding season for eight species: American Goldfinch, American Robin, Blue Jay, Common Grackle, European starling, House Finch, Purple Finch, and White-throated Sparrow (Figure 1E, 2D; Appendix C). The effect was non-significant for the other species: positive in two cases, and negative in two (Appendix C).

Discussion

Here we investigated whether natural and anthropogenic conditions affect the propensity of individuals to remain resident in twelve partially migratory North American bird species, using a large scale ringing dataset (USGS Bird Banding Laboratory, 2016).

In all species studied, we found substantial individual variability, with the same conditions under which some individuals migrate apparently tolerated by others year-round. This may reflect true individual variation in migratory strategies within populations. For example, individual differences in propensity to migrate have been related to factors such as sex [START_REF] Dhondt | Population Dynamics and Partial Migration of the European Robin (Erithacus rubecula) in Different Habitats[END_REF][START_REF] Perez | When males are more inclined to stay at home: insights into the partial migration of a pelagic seabird provided by geolocators and isotopes[END_REF], dominance status [START_REF] Ketterson | Seasonal, Annual, and Geographic Variation in Sex Ratio of Wintering Populations of Dark-Eyed Juncos (Junco hyemalis)[END_REF], personality (reaction to a novel object; [START_REF] Nilsson | Migratory and resident blue tits Cyanistes caeruleus differ in their reaction to a novel object[END_REF], body size [START_REF] Belthoff | Partial migration and differential winter distribution of House Finches in the eastern United States[END_REF] and physiology (e.g. basal metabolic rates and cost of thermoregulation; [START_REF] Nilsson | Basal metabolic rate and energetic cost of thermoregulation among migratory and resident blue tits[END_REF]. Individual strategies can also change over time, for example in response to changing environmental conditions [START_REF] Shaw | To breed or not to breed: a model of partial migration[END_REF] or age, an effect that we could not disentangle from individual variability, since we only had one pair of observations per individual.

It is also possible that this individual variability arose from limitations in our data, for instance in the environmental axes we considered and their proxies. For example, NDVI as a general measure of resources does not necessarily capture the specific resources needed for each species. And human density is a crude proxy for anthropogenic effects (e.g. agricultural areas often have low population densities and yet can provide important resources, [START_REF] Foley | Solutions for a cultivated planet[END_REF]; managed green spaces are often more productive than the surrounding wildlands, Imhoff et al. 2009).

Apparent individual variability in the propensity to migrate for apparently similar conditions can also arise from the temporal mismatch between the ringing data (see Appendix B for a distribution of records through time) and the explanatory variables (e.g. temperature averaged over 1970-2000, population density in 1970). Furthermore, data limitations meant we were unable to integrate withinseason mobility (e.g. [START_REF] Thorup | Resource tracking within and across continents in long-distance bird migrants[END_REF], yearly variation in migratory propensity (e.g. species known for their irruptive migrations: Evening Grosbeak, , [START_REF] Bock | Synchronous Eruptions of Boreal Seed-Eating Birds[END_REF]Pine Siskin, Alsop 2002), and the possibility that among our study species there may have been changes in migratory propensity over time (e.g. migratory populations becoming increasingly resident, [START_REF] Adriaensen | Ringing recoveries and the increase in numbers of European Great Crested Grebes Podiceps cristatus[END_REF][START_REF] Brown | Band recoveries reveal alternative migration strategies in American Robins[END_REF] or the opposite, Berthold 1999). Despite these limitations, which likely added noise to our data, our results support the predictions that local environmental conditions as we measured them affect the migratory decisions of individuals, in agreement of previous studies at the species level [START_REF] Somveille | Why do birds migrate? A macroecological perspective[END_REF]. In particular, and for all but one species, our results strongly support the hypothesis that winter harshness drives individuals to migrate elsewhere from their breeding locations (Fig 1B ,2A). For most species, we also found support for the hypothesis that high local natural surpluses in summer increase the propensity of individuals to remain as residents in their wintering locations, even if we only found a significant effect in five cases (Fig 1D ,2C).

For nine species, results also support the hypothesis that urbanization (measured through human population density) affects the propensity of individual birds to migrate (Fig. 1C,1E,2B,2D; Appendices B and C): in four species by increasing the probability that individuals stay over winter in their breeding locations, possibly by buffering them against winter harshness; in eight species by increasing the likelihood that individuals remain during the breeding season in their wintering grounds, possibly by increasing local resources. Our analysis thus adds to the existing evidence that urbanization is contributing to sedentarize at least some migratory species, showing that this effect can happen in either the wintering or the breeding portions of species' ranges.

It is not immediately clear why some species appeared in our results to be more responsive to urbanization than others. This could be in part due to limited statistical power, as the species for which we found no effect (Brown-headed Cowbird, Pine Siskin and Red-winged Blackbird) are among the ones with the smallest number of records (respectively 131, 104 and 160; but we found an effectin the winter onlyfor the White-throated Sparrow and the House Finch, with respectively 69 and 123 records). Species with irruptive migrations (i.e. which move in irregular patterns in the winter tracking blooms in resources), may be less responsive to human density. Indeed, we found that the two irruptive migrants in our sample, Pine Siskins and Evening Grosbeaks, showed no response to human density in their propensity to remain at their winter location (Fig. 2B). However, Evening Grosbeaks seem to respond to human density in their propensity to remain at their breeding location (Fig. 2B).

Another explanation may come from species' level of association with anthropogenic resources and infrastructures. Although a review by [START_REF] Archer | A systematic review of forest bird occurrence in North American forest fragments and the built environment[END_REF] found strong support for most species in our sample being synanthropic (medium support for the American Goldfinch, low support for the American Robin, and a lack of data about the White-throated Sparrow), there are known differences between species. For example, North-American populations of European Starlings are highly associated with cities [START_REF] Alsop | Smithsonian Handbooks: Birds of New England[END_REF], and we found a significant positive effect of human density on the propensity to remain as resident in both seasons (Fig. 1C, 1E, 2B, 2D; Appendices B and C). Similarly, American Goldfinches are common in suburbs, parks and backyards, and tend to visit garden feeders [START_REF] Alsop | Smithsonian Handbooks: Birds of New England[END_REF], and we found for this species a higher propensity to be resident (in both seasons) in densely human-populated areas (Fig. 2B,2D). And House Finches in the east of North American are common in human-created habitats [START_REF] Hill | A red bird in a brown bag: the function and evolution of colorful plumage in the House Finch[END_REF] and were found to have a higher propensity to remain in their wintering location over the summer in densely human-populated areas (Fig. 2D). In contrast, we found no effect of human density for Red-winged Blackbirds, whose main habitats (marshes and agricultural fields) are outside densely populated areas [START_REF] Alsop | Smithsonian Handbooks: Birds of New England[END_REF].

A higher propensity to remain year-round in urban areas does not obviate the fact that urbanization itself often has substantial impacts on bird diversity, resulting in poorer communities [START_REF] Lee | Breeding bird species richness in Taiwan: distribution on gradients of elevation, primary productivity and urbanization[END_REF][START_REF] Biamonte | Effect of urbanization on the avifauna in a tropical metropolitan area[END_REF]. In particular, a meta-analysis showed a lower passerine fledging success in urban areas [START_REF] Chamberlain | Avian productivity in urban landscapes: a review and meta-analysis[END_REF]. Nonetheless, our results suggest that, at least for some bird species that can cope well with anthropogenic habitats, urban areas can provide conditions that are beneficial enough to outweigh the costs of migration. However, it is possible that the benefits of urban areas only stand up to a certain level of urbanization, in the same way as Tratalos et al. ( 2017) found a hump-shaped relationship between bird species richness or abundances and household densities in Britain. Sample sizes for some of our study species were too small to allow us for the inclusion of nonmonotone effects in our models, so we were unable to test for this additional hypothesis.

Additionally, it is unclear whether individual choices to remain in urban areas rather than migrating translate into positive effects at the population level. Indeed, urban areas may be acting as ecological trapsi.e. when organisms choose poor-quality habitats above better alternatives (Gilroy and Sutherland 2006). Urban resources may appear more attractive yet result in worse outcomes; for example, great tits have been shown to prefer larger nesting cavities in urban areas despite this leading to lower fledging success [START_REF] Demeyrier | Experimental demonstration of an ecological trap for a wild bird in a human-transformed environment[END_REF]. Concerns have also been raised about the nutritional value of anthropogenic food [START_REF] Jones | Feeding birds in our towns and cities: a global research opportunity[END_REF], the increased risk of disease spread around anthropogenic food sources (Robb et al. 2008;[START_REF] Jones | Feeding birds in our towns and cities: a global research opportunity[END_REF]. Previous studies suggest complex effects: for example, for garden feeders, a review by Robb et al. (2008a) found mostly positive effects on the breeding performance of supplementary-fed birds (e.g. Robb et al. 2008b), but [START_REF] Plummer | Winter food provisioning reduces future breeding performance in a wild bird[END_REF] found a negative effect of feeders on the breeding performance of blue tits. To estimate whether remaining as resident in urban areas is actually beneficial would require data on population trends for migratory versus non-migratory individuals, or monitoring their breeding performance, which we did not have access to in this study.

Our results also indicate that, for some species, the effects of urbanization on the propensity to migrate are of similar orders of magnitude as those of variation in natural environmental conditions (as measured by the standardised coefficients of the regressions; Fig 2). This suggests that recent human activities are changing environmental conditions at scales comparable to those naturally experienced by species over longer-term evolutionary-ecological processes. Based on our results, we predict that the combination of climate change (leading to an increase in average temperatures) and urbanization will converge to decrease the propensity of individuals of many migratory bird species across their current range in temperate areas. But for many of those, it will also lead to their expansion into new higher-latitude breeding grounds. Furthermore, the effects of climate change are spatially and temporally complex, as even though average temperature increase, some areas may experience a higher frequency of extremely cold winters. Climate change can also affect the spatial and temporal distribution of local resources, in particular through the interactions between temperature and precipitation on local primary productivity, which may increase in some areas and decline in others. Finally, humans are affecting the distribution of local resources through land use change, including not only urbanization but also other habitat changes at continental scales [START_REF] Foley | Global consequences of land use[END_REF]. Overall, our results add to the body of evidence that bird migration is not a fixed behaviour in bird populations, and that it is already being shaped by human activities. Numbers indicate sample sizes: number of residents / total number of birds. B) to E) Environmental conditions experienced by individuals, represented as density curves (area under the curve equals one). Each graph indicates the frequency distribution of two sets of individuals: residents (in orange) and migrants (blue for winter migrants; red for summer migrants). B) and C) correspond to environmental conditions in the breeding locations; D) and E) to conditions in the wintering locations. Hence, for example, for the American Robin: in B a strong negative deviation of the red curve (migrants) in relation to the orange curve (residents) indicates that individual birds that migrate away from their breeding locations (i.e. summer migrants, thus in red) tend to be found in locations that experience lower temperatures in the winter than the those occupied by individuals that remain year round (residents, in orange). This is then reflected in Figure 2A by a significant positive effect of the winter temperature at the breeding location on the propensity of individuals to migrate.

Figure 1 :

 1 Figure 1: Geographic and environmental distribution of bird records, for each of the twelve species

Figure 2 :

 2 Figure 2: Migratory response of individuals to environmental conditions, analysed through binomial Generalised Linear Models (GLMs) modelling the propensity of individuals to remain as resident as function of local conditions. Values correspond to estimated coefficients of the GLMs for each species with 95% confidence intervals, with positive values indicating a higher propensity to remain as resident, and negative values a higher propensity to migrate. Transparency: significance after Bonferroni correction, with coefficients significantly different from 0 in dark. A) Effects of winter temperatures at the breeding locations. B) Effects of urbanisation at the breeding locations. C) Effects of breeding resource surplus at the wintering locations. D) Effects of urbanisation at the wintering locations.

Supplementary materials

Appendix A: Temporal distribution of ringing / re-sighting events for each species. House finch records correspond to only the Eastern population (East of 95˚W). Purple: migrants; orange: residents.

Appendix B: Results of the model predicting the probability of individuals remaining year-round at their breeding locations. The model is a binomial GLM (logit link) with the resident (1) versus migrant (0) status of each individual as a response variable and local environmental conditions (winter temperature and human population density) as predictors. The table presents estimated coefficients and p-values for the two predictor variables, with p-values in bold whenever the effect is significant after a sequential Bonferroni correction, and percentage of deviance explained by the model with the two variables.

Species

Erratum

In this analysis, each individual bird was classified into either 'resident' or 'migrant' based on thresholds in the distance between its breeding and non-breeding locations. In the published article, we defined these thresholds as: residents <20km, migrants >100km. However, a subsequently detected coding error means that the thresholds that were in practice used in the analyses were: residents <51.8km, migrants >259.0km.

No main results or conclusions were affected by these changes.


Marina Selezneva

Naïm Naouar

Yvan Denis

Larissa Gorbatikh

Peter Hine

Stepan V Lomov
email: stepan.lomov@kuleuven.be

Yentl Swolfs

Ignaas Verpoest

Philippe Boisse

Naim Naouar

Identification and validation of a hyperelastic

Keywords: self-reinfroced polypropylene, draping, simulations, material properties, wrinkles

come

Introduction

Polypropylene, which is a commodity polymer, can be used to produce tapes that have high molecular orientation and higher mechanical properties than the standard non-oriented polypropylene. These tapes when combined with polypropylene matrix, which is non-oriented, constitute a self-reinforced polymeric composite (SRPP). The composite is produced by heating up the tapes, which causes their outer sheath to melt while the inner core maintains its orientation. During the cool-down, the molten polymer consolidates and forms the "matrix" component of the SRPP. As these composites are reheated for forming, the oriented molecules can shrink and partially lose their orientation [START_REF] Alcock | Single polymer composites based on polypropylene: processing and properties[END_REF]. Hence, SRPP preforms need to be constrained during forming to avoid excessive shrinkage. Cabrera et al. [START_REF] Cabrera | Non-isothermal stamp forming of continuous tape reinforced all-polypropylene[END_REF] investigated deformation modes of woven SRPP during stamp-forming. They found that behaviour of SRPP is similar to that of the other textile-based composites and the main deformation mechanism is intra-ply shear [START_REF] Cabrera | Non-isothermal stamp forming of continuous tape reinforced all-polypropylene[END_REF]. However, SRPP was also shown to undergo additional tape drawing, which is a mechanism particular to self-reinforced polymeric composites and may be beneficial when forming complex geometries. The effect of the clamping conditions on formability of SRPP composites was investigated experimentally by Selezneva et al. [START_REF] Selezneva | Formability Study of Hybrid Carbon Fibre/Self-Reinforced Polypropylene Composites[END_REF]. Three clamping options were considered: full edge, spring and corner. These conditions triggered different degrees of stretching, shearing and draping of the preforms. As the result, the formed hemispheres exhibited different wrinkling patterns.

Draping modelling is an indispensable step of a thermoformed composite part design. There are two main factors, related to the draping, which affect mechanical behaviour of the final consolidated part.

International Journal of Material Forming (2020) DOI: 10.1007/s12289-020-01542-3

First, the "regular" change of the local composite properties because of the draping. The draping-induced reinforcement deformation, primarily shear, changes local geometry of the fibres. The basic mechanical properties of the composite, such as stiffness and strength, therefore can strongly change due to the draping. This is a topic that's been studied a lot, especially in recent years [START_REF] Aridhi | Textile composite structural analysis taking into account the forming process[END_REF][START_REF] Boisse | Composites forming[END_REF][START_REF] Boisse | Composite Reinforcements for Optimum Performance[END_REF][START_REF] Gereke | Experimental and computational composite textile reinforcement forming: A review[END_REF][START_REF] Lomov | Micro-macro structural analysis of textile composite parts[END_REF][START_REF] Tavana | Meso/macro-scale finite element model for forming process of woven fabric reinforcements[END_REF][START_REF] Van Den Broucke | Determination of the mechanical properties of textile-reinforced composites taking into account textile forming parameters[END_REF][START_REF] Davidson | The effects of defects on the compressive response of thick carbon composites: an experimental and computational study[END_REF][START_REF] Han | Draping simulation of carbon/epoxy plain weave fabrics with nonorthogonal constitutive model and material behavior analysis of the cured structure[END_REF]. These analyses are included in integrated composite part design software (see, for example [START_REF] De Luca | Numerical simulation of reinforcements forming: the missing link for the improvement of composite parts virtual prototyping[END_REF][START_REF] Bayraktar | Forming and performance analysis of a 3D-woven composite curved beam using meso-scale FEA[END_REF][START_REF] Kärger | Forming optimisation embedded in a CAE chain to assess and enhance the structural performance of composite components[END_REF]). Once the constitutive model of reinforcement deformation during draping is identified, the local fibres geometry can be modelled quite accurately with the state-of-the-art draping modelling software, and the local mechanical properties can be calculated [START_REF] Lomov | Meso-level textile composites simulations: open data exchange and scripting[END_REF][START_REF] Verpoest | Virtual textile composites software Wisetex: integration with micromechanical, permeability and structural analysis[END_REF].

Second, defects of the draping, the most prominent of which is wrinkling, can lead to drastic deterioration of the load-carrying ability of the consolidated part and can cause part rejection. Some studies have analysed the influence of the presence of wrinkles, created during manufacturing, on the mechanical performance of the composite [START_REF] Leong | Fatigue failure of sandwich beams with face sheet wrinkle defects[END_REF][START_REF] Bloom | Damage progression and defect sensitivity: an experimental study of representative wrinkles in tension[END_REF][START_REF] Mukhopadhyay | Compressive failure of laminates containing an embedded wrinkle; experimental and numerical study[END_REF][START_REF] Bender | Investigation of the effect of wrinkle features on wind turbine blade sub-structure strength[END_REF]. A recent review [START_REF] Boisse | Bending and wrinkling of composite fiber preforms and prepregs. A review and new developments in the draping simulations[END_REF] demonstrates that the presence, pattern and intensity of wrinkling can be accurately modelled once the constitutive model of the reinforcement is identified. The part design can be then done "as manufactured", accounting for the draping process induced defects.

Mechanical behaviour of textile reinforcements during draping is an important subject in composite materials science, with hundreds of papers published. Among these papers not that many are dedicated to the identification of constitutive models of thermoplastic pre-impregnated textile sheets at the processing temperature. Since the benchmark exercise [START_REF] Cao | Characterization of mechanical behavior of woven fabrics: Experimental methods and benchmark results[END_REF] there was a considerable amount of research done on glass fabric reinforced polypropylene for quantification of the constitutive models (for example, [START_REF] Kuhtz | Characterising the Thermoforming Behaviour of Glass Fibre Textile Reinforced Thermoplastic Composite Materials[END_REF][START_REF] Machado | Analysis of the thermomechanical shear behaviour of woven-reinforced thermoplastic-matrix composites during forming[END_REF][START_REF] Yin | Forming of thermoplastic plain woven carbon composites: An experimental investigation[END_REF]). The forming studies of SRPP [START_REF] Abdiwi | Measuring the Shear-Tension Coupling of Engineering Fabrics[END_REF][START_REF] Alcock | Direct Forming of All-Polypropylene Composites Products from Fabrics made of Co-Extruded Tapes[END_REF][START_REF] Kalyanasundaram | Effect of process parameters during forming of self reinforced -PP based Fiber Metal Laminate[END_REF][START_REF] Kalyanasundaram | Effect of temperature and blank holder force on non-isothermal stamp forming of a self-reinforced composite[END_REF][START_REF] Kim | Optimum consolidation of self-reinforced polypropylene composite and its time-dependent deformation behavior[END_REF][START_REF] Zanjani | The Effect of Fiber Orientation on the Formability and Failure Behavior of a Woven Self-Reinforced Composite[END_REF] are focused on determining process windows for forming of SRPP in the presence of shrinkage rather than on identification of the constitutive behaviour at processing temperature.

The present work focuses on the development and validation of an FE model for thermoforming of SRPP. The modelling approach is based on the hyperelastic law developed by Charmetant et al. [START_REF] Charmetant | Hyperelastic model for large deformation analyses of 3D interlock[END_REF]. The constitutive model considers the main deformation modes that occur during the forming of textile composites, such as stretching in the warp and weft directions and in-plane shear. Material behaviour associated with each deformation mode is defined in terms of a physical invariant and a strain energy density function. In the study, the strain energy functions were obtained by fitting experimental data from tensile and picture-frame shear tests, which were conducted at forming temperatures. To validate the FE model, results were compared against the experimental data for the spring-clamping condition. The states or geometries of the SRPP preform at different stages of the forming process were captured by moulding shallow domes with depths ranging from 1-4 cm by using a hemispherical mould with a 4.5 cm radius. To the authors' knowledge, this is the first work where the parameters of the Boisse-Charmetant formulation are identified for SRPP and are validated in comparison with experiment, hence can be used by others for virtual investigation of the SRPP forming.

In this paper, first the material and the experimental forming studies are discussed. Next, the description of the constitutive model and the development of the FE simulation are presented. Finally, results and conclusions are given.

Experimental

Material and Preform Consolidation

The precursor material to SRPP is a balanced 2/2 twill PP tape fabric provided by Propex Fabrics GmbH, as shown in Figure 1a. The preform panels were prepared by stacking in an alternating manner twelve PP fabrics and eleven 20 µm thick PP films, as shown in Figure 1b. Panels were consolidated by hot compaction at 188°C with a hold time of 5 min at 40 bar and an average cooling rate of 35°/min. This process was optimized for SRPP with tensile properties in mind [START_REF] Swolfs | The influence of process parameters on the properties of hot compacted self-reinforced polypropylene composites[END_REF]. The produced panels were 22 cm by 22 cm and 2 mm thick.

Mechanical Testing

To measure the tensile and in-plane shear behaviour at forming temperatures, tests were performed on 2 mm thick pre-consolidated panels in a convection oven at 170°C.

Picture frame tests, for in-plane shear properties, were performed using the picture frame setup developed at KU Leuven, and mounted on an Instron 5985 machine with a 100 kN load cell. The geometry of the picture frame fixture is shown in Figure 2a, the picture frame is installed vertically on a tensile frame. Additional information regarding the experimental set up and the data processing associated with this test can be found in [START_REF] Lomov | Picture Frame Test of Woven Composite Reinforcements with a Full-Field Strain Registration[END_REF]. The picture frame setup was placed in a convection oven, and had to be placed at 45° to fit within the oven. Unfortunately, this prevented the use of digital image correlation. The tests were performed at a displacement rate of 20 mm/min, which corresponded to a shear angle rate of approximately 50°/min (decreasing from 54°/min at the beginning of the test to 45°/min at the maximum machine stroke, as calculated with formulae describing the kinematics of the frame in [START_REF] Lomov | Picture Frame Test of Woven Composite Reinforcements with a Full-Field Strain Registration[END_REF]). Six calibration tests with an empty frame were performed to account for friction in the setup. The average of these force-displacement diagrams was used to subtract the frictional force from the test results. After this correction, the shear angle and stress were calculated using the formulas described in [START_REF] Lomov | Picture Frame Test of Woven Composite Reinforcements with a Full-Field Strain Registration[END_REF]. Tensile samples had a width of 20 mm and a gauge length of 80 mm, which was limited by the dimensions of the oven. Tests were carried out at 20 mm/min, which is equivalent to a nominal strain rate of 25%/min. Tests were interrupted when the samples began to pull out of the grips. Tensile tests were affected by the transverse shrinkage of the SRPP. Therefore, for future work it is recommended to test self-reinforced polymer composites using a biaxial test. This test would restrain the sample better and should minimize shrinkage. Furthermore, the induced deformation state in a biaxial test is also closer to the real deformation state during thermoforming.

As was discussed by Liang et al. [START_REF] Liang | Analysis of thermoplastic prepreg bending stiffness during[END_REF] and Boisse et al. [START_REF] Boisse | Bending and wrinkling of composite fiber preforms and prepregs. A review and new developments in the draping simulations[END_REF], bending stiffness is another important property to be included in the simulation, as it dictates the number and size of wrinkles. It has been shown that lower stiffness leads to finer wrinkles [START_REF] Liang | Analysis of thermoplastic prepreg bending stiffness during[END_REF]. In contrast to what is done for continuous materials such as metals, the bending International Journal of Material Forming (2020) DOI: 10.1007/s12289-020-01542-3 rigidity of textiles cannot be calculated directly from the tensile modulus [START_REF] Boisse | Bending and wrinkling of composite fiber preforms and prepregs. A review and new developments in the draping simulations[END_REF]. Bending rigidity of textiles can be measured using the Peirce cantilever test [START_REF] Peirce | The "handle" of cloth as a measurable quantity[END_REF]. It is based on the cantilever bending of a textile specimen under its own weight, as shown in in Figure 2b. The specimen is progressively advanced until the free end makes contact with the inclined plane of the device (Figure 2b). Assuming a linear relation between the bending moment M and the curvature , the specific flexural rigidity S is determined from the overhanging length l corresponding to an inclination of 41.5°:

M = G S G w 3 8

where G is the bending stiffness of the specimen, w is the weight per unit length of the specimen.

Tests on fabrics can be conducted at ambient temperatures, but composites should be tested at forming temperature to account for the matrix softening. However, SRPP composite and PP fabric would shrink at forming temperatures. Hence, PP fabric was tested at room temperature and the results were treated as the upper limit of the bending rigidity. In further refinement of the model, the viscosity of the molten polypropylene can be accounted for in addition to the bending resistance of the fabric [START_REF] Harrison | Rate dependent modelling of the forming behaviour of viscous textile composites[END_REF].

The loading speed during the forming was 20 mm/s, hence the hemisphere was formed in 2-3 s. This gives an estimation of the maximum local shear rate of the fabric of 5-10°/s (referring to the maximum shear angle of 20°, see Error! Reference source not found.), with the minimum shear rate close to zero. The shear rate in picture frame tests was 0.7-0.9°/s, as stated in the previous section. Therefore in most of the hemisphere surface the shear rates in the forming and in the picture frame tests are of the same order of magnitude. The one decimal order of magnitude difference in the zones of maximum shear can indeed affect the results; this effect should be estimated in future work.

Thermoforming

The thermoforming tests were performed at the facilities of the SLC Lab (Leuven). A Fontijne press with integrated infrared heaters were used for forming. Two types of boundary conditions were considered: springsupport and full-edge clamping (actually the clamping was not ideal and certain blank sliding occurred). These conditions triggered different degrees of stretching, shearing and draping of the preforms. Error! Reference source not found. shows a schematic of the spring-supported and fully-clamped set up. The springs had a stiffness of 260 N/m. Preforms used for forming were 2 mm thick and 22 x 22 cm in size. Preforms were heated using infrared heaters. Temperature was monitored using embedded thermocouples (in the middle and near the surface) and the pyro-sensors incorporated into the heaters. Material were heated to a surface temperature of 170-175°C. Once the set temperature was reached, it was transferred to the press using a rail system and stamped. The mould was kept at 100°C, and the panels were demoulded once their temperature cooled down and stayed at 100°C for 2-3 minutes. Throughout the study, a steel hemispherical mould with a 4.5 cm radius was used. To capture how the material behaves during forming, shallow domes were formed by partially closing the mould to 1, 2, 3 and 4 cm depths. The material was then allowed to cool down and re-solidify while preserving this intermediate geometry.

International Journal of Material Forming (2020) DOI: 10.1007/s12289-020-01542-3

3 Constitutive model and parameters identification

Hyperelasticity

The hyperelastic constitutive model developed by Charmetant et al. [START_REF] Charmetant | Hyperelastic model for large deformation analyses of 3D interlock[END_REF] was used to model the highly nonlinear mechanical properties of SRPP. The required input parameters were fitted using the experimental data. In the hyperelastic framework, the Cauchy stress tensor (σ) can be calculated from the deformation gradient tensor (F) and the strain energy potential (w):

σ 2 J F • δw δC • F ( 1
)
where C is the Cauchy-Green deformation tensor and J is the Jacobian determinant of F. Material properties of the initially orthotropic material are defined along the three orthogonal directions: warp M 1 , weft M 2 and through-the-thickness M 3 . Hence, strain energy density function for a hyperelastic law can be expressed [START_REF] Quanshui | Tensor function representations as applied toformulating constitutive laws for clinotropic materials[END_REF][START_REF] Itskov | On the theory of fourth-order tensors and their applications in computational mechanics[END_REF] in terms of invariants: and where the mixed invariants correspond to the structural tensors:

I M • C • M 𝑖 1,3 I M • C • M 𝑖 1,3 I M • C • M 𝑖 1,3
A constitutive model for initially orthotropic materials can be derived on the basis of elementary deformation mode (e.g. in-plane elongation, shear) [START_REF] Charmetant | Hyperelastic model for large deformation analyses of 3D interlock[END_REF]. It can be assumed that the deformations modes are uncoupled and that the total strain energy density function can be written as the summation of n strain energy densities. Hence, the derivative of w can be expressed as:

δw δC δw δI δI δC ( 2
)
where I i is the strain invariant of the ith deformation mode. The expressions for the invariants for each deformation mode were derived by Charmetant et al. [START_REF] Charmetant | Hyperelastic model for large deformation analyses of 3D interlock[END_REF] based on physical observations. The main deformation modes influencing wrinkling of textiles are in-plane elongation, in-plane shear and transverse shear, as shown in Figure 4. The method for obtaining the strain energy function for these modes using the experimental data will be described hereafter.

In-plane elongation

For in-plane elongation along the warp and weft fibres, Charmetant et al. [START_REF] Charmetant | Hyperelastic model for large deformation analyses of 3D interlock[END_REF] proposed the following expressions for the elongation invariant (I elong ) and its derivative:

I ln I ε (3) δI δC 1 2I M
In case of uniaxial loading, I elong equals to the experimentally measured tensile strain (ε), and Eq. ( 1) for the Cauchy stress can be rewritten as:

International Journal of Material Forming (2020) DOI: 10.1007/s12289-020-01542-3

σ 2 J F • δw δC • F 2 J F • δI δC • F δw δI
The Cauchy stress component corresponding to the longitudinal tension is hence:

σ 1 I δw δI (4)
Experimental data is fitted by working with the curve for I elong vs. dw elong /dI elong calculated using Eq. ( 3) and ( 4), respectively. The fitted expression for dw elong /dI elong can then be integrated to obtain the coefficients of the strain energy density (w) function. The piecewise polynomial form of the energy function that was obtained is shown in Eq. ( 5), and the tensile stiffness coefficient K t0 and K t are given in Table 1. The correlation between the experimental and fitted data is shown in Error! Reference source not found.a.

w ⎩ ⎨ ⎧ 1 2 K I 1 2 K I 1 2 K I 0.261 0.261 K I 1 2 K I 𝐼 0 𝐼 0.261 𝐼 0.261 (5)
No experimental data was available for the compressive behaviour (K c ) of SRPP at elevated temperature. It was previously measured that at room temperature K c ≈ K t [START_REF] Mesquita | Flexural behaviour of corrugated panels of self-reinforced polypropylene hybridised with carbon fibre: An experimental and modelling study[END_REF]. However, at elevated temperature SRPP shrinks along the warp and weft axis, and hence is expected to have low compressive modulus. Charmetant et al. [START_REF] Charmetant | Hyperelastic model for large deformation analyses of 3D interlock[END_REF] assumed that K c of 3D interlock textile is equal to the tensile stiffness at low strains (<1.5%), which was significantly lower than stiffness at higher strains (38 vs. 816 MPa). On the other hand, Pazmino et al. [START_REF] Pazmino | Numerical modelling of forming of a non-crimp 3D orthogonal weave[END_REF] warned that a large mismatch between compressive stiffness and transverse shear rigidity can cause nonphysical wrinkling, as shown in Error! Reference source not found.b. To avoid this problem, some authors set K c = 0.01K t for non-crimp 3D textiles [START_REF] Pazmino | Numerical modelling of forming of a non-crimp 3D orthogonal weave[END_REF]. In the present study, two values of K c were used in the simulation (K c = K t0 and K c = 0.1K t0 ), shown in Error! Reference source not found.a. As will be discussed in the results section, value that gave the best correlation with experimental results was recommended for future use. International Journal of Material Forming (2020) DOI: 10.1007/s12289-020-01542-3

In-plane shear

For in-plane shear, Charmetant et al. [START_REF] Charmetant | Hyperelastic model for large deformation analyses of 3D interlock[END_REF] proposed the following expressions for the shear invariant (I shear ) and its derivative:

I I I I sin γ (6) δI δC 1 2 I I M ⊗ M M ⊗ M I 2I M I 2I M
where γ is the shear angle. The Cauchy stress component for in-plane shear (σ 12 ) is thus expressed as:

σ δw δI cos γ (7)
The experimentally measured curve is shown in Figure 6a. This curve is atypical since there is no region of fast increase in the stress values that would correspond with the locking angle. However, from the test video it was obvious that the specimen developed wrinkles during the picture frame test, as shown in Figure 6b. Hence, it was decided to fit the initial part of the curve, while assuming that the second half corresponds to the mixed mode shearing and wrinkling or bending behaviour. The resultant fitted curve is defined in two parts to ensure that the computed stress is negative for negative values of I shear , as shown in Eq. ( 8). To check the effect of this assumption, simulations were also run using the fully fitted experimental curve, and the resultant expression for w shear is given by Eq. ( 9). The units for I shear and w shear are mm/mm and MPa.

w K • I K • I K • I K • I 𝐼 0 𝐼 0 (

Transverse shear

Bending stiffness is an important parameter as it dictates the shape of wrinkles [START_REF] Liang | Analysis of thermoplastic prepreg bending stiffness during[END_REF]. However, the 3D solid elements that were used in this simulation do not have bending rigidity. Bending properties can, however, be imposed indirectly via transverse shear. The governing equations for transverse shear are defined in a similar way as for in-plane shear, hence transverse shear invariant (I trans ) and its derivative are:

I I I I si n γ , α 1,2 δI δC 1 2 I I M ⊗ M M ⊗ M I 2I M I 2I M , α 1,2 w K I (10)
Transverse shear behaviour at room temperature was measured indirectly by a Pierce cantilever bending test [START_REF] Pazmino | Numerical modelling of forming of a non-crimp 3D orthogonal weave[END_REF], as shown in Figure 2b. An FE model was created using the same geometry and boundary conditions, and transverse shear was adjusted until the same vertical displacement was obtained. In these calculations, one element in the model thickness was used, to have an isolated influence of the transversal shear on the bending resistance. The polynomial equation that was used to represent the transverse shear behaviour is shown in Eq. [START_REF] Van Den Broucke | Determination of the mechanical properties of textile-reinforced composites taking into account textile forming parameters[END_REF]. The same form was adopted by Pazmino et al. [START_REF] Pazmino | Numerical modelling of forming of a non-crimp 3D orthogonal weave[END_REF]. Based on the cantilever test conducted at room temperature, the transverse stiffness coefficient K trans is about 5 MPa. This value represents the upper limit of stiffness, since PP fibres will soften at high temperatures and lose their bending rigidity. The coefficient K trans was finally adjusted in the dome forming model to achieve the correct wrinkle geometry, as will be discussed in the results section.

Model limitations and heuristic assumptions

The use of the hyperelastic model, described above, involves a number of heuristic assumptions.:

1. An important assumption is that the material model parameters, identified in the tests described above and performed with a single layer of the reinforcement, can be used as parameters of 3D elements in the stack of 12 layers. This assumption is especially important for in-plane shear and bending (transverse shear). 2. The viscoelastic nature of the material resistance is not included in the model. The approach in [START_REF] Harrison | Rate dependent modelling of the forming behaviour of viscous textile composites[END_REF] can be a basis for future development of a coupled visco-hyperelastic model. Certain justification for use of International Journal of Material Forming (2020) DOI: 10.1007/s12289-020-01542-3

(hyper)elastic model is given by a recent work [START_REF] Ghafour | The importance of taking into account behavior irreversibilities when simulating the forming of textile composite reinforcements[END_REF], which has compared simulations of composite forming using an elastic model on the one hand and irreversible models on the other. It has been shown that in the cases of monotonous forming the results are close in the both cases. A rigorous verification of the validity of these assumptions and further model development will be a subject of future work.

Simulation

The finite element model was created in Abaqus 6.14. By taking advantage of the symmetry conditions, a quarter of the panel geometry was modelled. The position of the quarter panel with respect to the rigid mould as well as the imposed symmetry conditions are shown in Figure 7. The SRPP preform was modelled using 3D explicit solid elements. Use of 3D solid elements allows full inclusion of transverse shear effects to represent bending, which is of major importance for wrinkling prediction [START_REF] Boisse | Bending and wrinkling of composite fiber preforms and prepregs. A review and new developments in the draping simulations[END_REF]. A mesh size of 1 x 1 mm, and four elements through the thickness were used. Behaviour of SRPP was modelled using a user-defined hyperelastic material (VUMAT), which was developed by Charmetant et al. [START_REF] Charmetant | Hyperelastic model for large deformation analyses of 3D interlock[END_REF].

The mould was modelled as a rigid body. Two types of boundary conditions were considered. In the clamped configuration, the outer edges of the panel were fully fixed. In the case with spring clamping, boundary conditions in the model were simplified by omitting the springs. This simplification is based on the assumption that springs offer low resistance, and the preform can slide freely into the mould.

Results and discussion

The distribution of the local shear angles of the fabric for a hemispherical draping is well defined by the kinematics of the draping and is not much influenced by the details of the constitutive mechanical model. Therefore the focus in the discussion of the results will be on the wrinkling behaviour, which is strongly influenced by the details of the material behaviour and by the relation between different components of the material resistance to the deformation [START_REF] Boisse | Bending and wrinkling of composite fiber preforms and prepregs. A review and new developments in the draping simulations[END_REF].

Simulation and test results for the (presumably) fully clamped specimen are shown in Figure 8. The stiffness parameters used in this simulation are summarized in Table 1. Location of some of the wrinkles is indicated by the red circles. While minor wrinkles are present in both experimental and simulation results, their location is different. In the case of the experimental samples, the edges came loose from the clamps as indicated by the arrow in Figure 8a. Hence, it can be postulated that the loss of tension in the middle of the edge permitted wrinkle formation in the highly sheared region. For future tests, the clamping fixture needs to be improved to achieve better clamping. Even if the clamping in the present experiments is not ideal, the results show sensitivity of the qualitative wrinkling characteristics to the boundary conditions. Comparison between the experimental and simulation results for spring-clamping is shown in Figure 9 and Figure 10. The stiffness parameters used in this simulation are summarized in Table 1. Experimental results are slightly non-symmetrical because of the difficulties associated with the alignment of the panel when spring clamping was used. Both experiments and calculations (see Figure 9a and Figure 10a) show that there are no wrinkles at the early stages of deformation. Once the material is deformed by 2 cm (see Figure 9b and Figure 10b), wrinkles all around the circumference begin to appear. When 4 cm of deformation is attained (see Figure 9c and Figure 10c), clear folds can be observed at the four sides. Overall, good qualitative correlation between experimental and modelling results is observed. Effect of the shear and compressive stiffness on wrinkling was also analysed to assess the significance of the assumptions. Images depicting the in-plane shear angle distribution for the two clamping conditions are shown in Error! Reference source not found.. In case of spring-clamping, in-plane shear strains were low (I shear ≤ 0.16). In this range, the fitted curves essentially overlap, refer to Figure 6a, and hence no difference was expected. In case of full clamping, the maximum shear angle was high (≈ 20⁰) and in the range of the imposed locking angle. Nonetheless, no wrinkles were predicted in that region regardless of which fitted curve was used. This can be explained by the high elongation strains created by the boundary conditions. Overall, in-plane shear stiffness had minimal effect on the wrinkles observed in both clamping cases. In the future it would be interesting to measure the shear angle of the formed pieces using a photogrammetry technique to advance the present qualitative comparison to quantitative level. [START_REF] Pazmino | Numerical modelling of forming of a non-crimp 3D orthogonal weave[END_REF] when there is a large mismatch between these two properties. When compressive modulus is decreased (moving to the right in figure), the shape of the wrinkles becomes more realistic. However, if it is decreased too much, small wrinkles smoothen out and disappear. An alternative approach, would be to keep the same compressive modulus and to instead increase the shear modulus (moving down in the figure), while keeping in mind that K trans = 5 MPa is the upper limit.

Conclusions

The parameters of the Boisse-Charmetant hyperelastic material model are identified for SRPP and are validated in comparison with experiment.

The parameters of the hyperelastic law are based on the experimental data from tensile and picture frame tests at the forming temperature. Flexural or transverse shear stiffness was adjusted to capture the correct wrinkle shape.

To further improve the model, this property will need to be measured experimentally. Another improvement in the model can be the inclusion of the interaction of the preform plies. Overall, the deformation of the preform during the thermoforming process was captured properly by the model. Additional improvements can be achieved

FIGURE 1 .

 1 FIGURE 1. (a) 2/2 twill PP woven fabric and (b) material layup used to produce preforms.

FIGURE 2 .

 2 FIGURE 2. (a) KU Leuven picture frame fixture: 1,2 -hinges, 3 -groove, 4 -lip, 5 -plate with screws (adapted from [23]), and (b) Peirce cantilever test fixture with clamped PP textile

FIGURE 3 .

 3 FIGURE 3. Experimental set-up with (a) spring clamping and (b) full-edge clamping

  w w I , I , I , I , I , I , I , I I I I where I 1 , I 2 and I 3 are the invariants of C:

FIGURE 4 .

 4 FIGURE 4. Main deformation modes of textiles during forming (adapted from[START_REF] Charmetant | Hyperelastic model for large deformation analyses of 3D interlock[END_REF]).

FIGURE 5 .

 5 FIGURE 5. (a) experimental and fitted tensile and compressive data, and (b) non-physical wrinkling (adapted from [42]).

8 ) 3 FIGURE 6 .

 836 FIGURE 6. (a) Experimental and fitted in-plane shear data and (b) wrinkled picture frame specimen at the end of the test.To check the effect of this assumption, simulations were also run using the fully fitted experimental curve, and the resultant expression for w shear is given by Eq. (9). The units for I shear and w shear are mm/mm and MPa. w 74 • I 192 • I 193 • I 91.4 • I 18.5 • I 0.65 • I 74 • I 192 • I 193 • I 91.4 • I 18.5 • I 0.65 • I

  FIGURE 6. (a) Experimental and fitted in-plane shear data and (b) wrinkled picture frame specimen at the end of the test.To check the effect of this assumption, simulations were also run using the fully fitted experimental curve, and the resultant expression for w shear is given by Eq. (9). The units for I shear and w shear are mm/mm and MPa. w 74 • I 192 • I 193 • I 91.4 • I 18.5 • I 0.65 • I 74 • I 192 • I 193 • I 91.4 • I 18.5 • I 0.65 • I 𝐼 0 𝐼 0 (9)

FIGURE 7 .

 7 FIGURE 7. Simulation geometry and boundary conditions.

FIGURE 8 .

 8 FIGURE 8. (a) Experimental and (b) simulation results of fully clamped preforms

FIGURE 9 .

 9 FIGURE 9. Experimental results of the intermediate forming stages of SRPP, spring clamping: (a) 1 cm high dome, (b) 2 cm high dome and (c) 4 cm high dome.

FIGURE 10 .

 10 FIGURE 10. Simulation results of the intermediate forming stages of SRPP, spring clamping: (a) 1 cm high dome, (b) 2 cm high dome and (c) 4 cm high dome.

FIGURE 11 .Figure 12 .

 1112 FIGURE 11. In-plane shear angle in simulations with (a) spring-clamping and (b) edge clamping. The effect of the compressive and transverse shear stiffness on the simulation results is shown in

Figure 12

 12 Figure 12. Examination of the left-top image (K trans = 0.1 and K c =89 MPa) reveals presence of non-physical wrinkles that were anticipated by Pazmino et al.[START_REF] Pazmino | Numerical modelling of forming of a non-crimp 3D orthogonal weave[END_REF] when there is a large mismatch between these two properties. When compressive modulus is decreased (moving to the right in figure), the shape of the wrinkles becomes more realistic. However, if it is decreased too much, small wrinkles smoothen out and disappear. An alternative approach, would be to keep the same compressive modulus and to instead increase the shear modulus (moving down in the figure), while keeping in mind that K trans = 5 MPa is the upper limit. The optimal combination of K c

FIGURE 12 .

 12 FIGURE 12. Effect of compressive and transverse shear modulus on wrinkle development in spring-clamped SRPP panel. All units are in MPa. The selected combination of K c and K trans is shown in bold.

TABLE 1 .

 1 Selected stiffness coefficients, MPa

	K t0	K t	K c =0.1K t0	K shear1	K shear2	K trans
	89	70	8.9	-0.18	16.875	0.1

DOI: 10.1007/s12289-020-01542-3 by incorporating the effect of the spring-clamping into the model. It would be interesting to further extend this model to hybrid composites, which have a broader spectrum of applications.

The comparison between the simulated and experimental wrinkling pattern shows that the identified constitutive model for SRPP draping represents the behaviour of the deformed sheet adequately and can be used for analysis of the influence of the process design choices (for example, the blank holder organisation) on the quality of the final part. The results of the local reinforcement deformation can be further transferred to the stress analysis of the consolidated part and will influence the stress-strain allowances based design.

Acknowledgements

MS is grateful to the ESAFORM mobility grant, which enabled the collaboration between INSA Lyon and KU Leuven to take place. The authors would also like to thank the SLC lab for their help with the forming trials. MS, LG and YSW gratefully acknowledge SIM (Strategic Initiative Materials in Flanders) and VLAIO (Flemish government agency for Innovation and Entrepreneurship) for their support of the ICON HYTHEC project. YSW also extends his gratitude to the FWO Flanders for his postdoctoral fellowship. The authors acknowledge the support provided by Toray Industries (SVL is the current holder of the Toray Chair for Composite Materials at KU Leuven).

Conflict of Interest Statement

The authors declare absence of any conflict of interests.