Beatrice Gleize

Dr Manon Hiolle

Dr Nathalie Meunier

Dr Bruno Pereira

Ruddy Richard

Dr Isabelle Savary-Auzeloux

Dr Caroline Buffière

Marie-Agnès Peyron

Dr Charlotte Halimi

Dr Catherine Caris-Veyrat

Béatrice Gleize

Françoise Nau

Emmanuelle Reboul
email: emmanuelle.reboul@univ-amu.fr

Dr B Gleize

Prof. F. Nau Stlo

Prof R Richard

Food Structure Modulates the Bioavailability of Triglycerides and Vitamin D, and Partly That of Lutein: A Randomized Trial with a Crossover Design in Adults

Keywords: carotenoids, clinical trial, cholecalciferol, chylomicrons, bioavailability

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Fat-soluble vitamins (i.e., vitamin A, D, E, and K) and bioactive phytochemicals (i.e., carotenoids) play key roles in both development and healthy aging. Although less frequent than in developing countries, fat-soluble vitamin (principally vitamin D and E) sub-deficiencies and/ or insufficient intakes have also been described in Western countries. [1,2] The link between these inadequacies and the increasing incidence of chronic diseases such as neurodegenerative, bone, or cardiovascular pathologies is now extensively investigated because many reports strongly suggest causative relationships. [3,4] The most frequent inadequacy concerns vitamin D 3 (cholecalciferol). [3] This vitamin is involved in many functions including the optimal functioning of the musculoskeletal system. [5] In most people, sun exposure is not sufficient to reach satisfactory vitamin D levels. [3] A complementary dietary intake is thus necessary. Optimizing the bioavailability of vitamin D 3 dietary intake would be of great interest. Another bioactive phytochemical of interest is the carotenoid lutein. This xanthophyll is not a "micronutrient" per se because its essentiality has not been proven in humans yet. Nevertheless, lutein was shown to selectively accumulate in the retina [6,7] where its presence has been negatively associated with the incidence of age-related macular degeneration, [8,9] which is a major cause of visual impairment in elderly. Lutein beneficial effect is likely due to its ability to protect photoreceptors from blue light. [10][11][12] Recent randomized, placebo-controlled clinical trials have confirmed these data as xanthophyll supplementation improves patient visual function and decreases the risk of progression to late age-related macular degeneration. [13] Optimizing dietary lutein intake would thus be of great interest as well.

Over the past decades, many studies explored the factors modulating fat-soluble micronutrient absorption (see for review [14,15] ). Indeed, even if these molecules are provided in sufficient amounts by the diet, their beneficial effects depend on their bioavailability, which can be affected by interactions with the matrix in which they are embedded.

The specific effect of food structure on the digestion of macronutrients such as proteins, [16] lipids, [17] or carbohydrates [18] has previously been explored, but studies on the effect of food structure on fat-soluble micronutrients remain scarce. In particular, there is no clear clinical data on fat-soluble vitamin release kinetics into the bloodstream after consuming foods with similar composition but different structures.

To fill this gap, we comprehensively compared both vitamin D 3 (as a model of fat-soluble vitamin) and lutein (as a model of carotenoid) postprandial responses in healthy subjects, after the consumption of four model foods presenting similar nutrient composition but different structures, that is, custard, pudding, sponge cake, and biscuit.

Experimental Section

Subject Number and Characteristics

The main objective of the study was to highlight the link between food structure and fat-soluble micronutrient postprandial responses. The primary outcomes were thus lutein and vitamin D 3 postprandial concentrations in chylomicrons. Previous works reported that 10-20 subjects would be relevant to perform kinetics studies having a satisfactory statistical power greater than 80%. [19,20] Twelve healthy, non-obese, non-smoking men were thus recruited for the study by the Unité d'Exploration Nutritionnelle (Centre de Recherche en Nutrition Humaine Auvergne, Clermont-Ferrand, France). Only men were enrolled to decrease interindividual variability in the observed responses. Subjects had no history of chronic disease, hyperlipidemia, or hyperglycemia and they were not taking any medication known to affect vitamin D 3 or lipid metabolism 3 months before the study or during the study period. They had to have at least 28 teeth, stable Angle Class I occlusion, no jaw pain, no current dental treatment, and no orthodontic procedure in the 3 years preceding the study. Subjects' baseline characteristics are reported in Table 1. Plasma glucose was evaluated at using hexokinaseglucose-6-phosphate dehydrogenase method. Hemoglobin was assayed by a quantitative and qualitative automated method (spectrometry, impedancemetry, cytofluorimetry); completed, if necessary, by cytological analysis and microscopic control. 25-hydroxyvitamin D status was evaluated with a chemilumines- cence immunological assay. Subjects were asked to avoid foods rich in the studied micronutrients for 48-h meals preceding each kinetics (see Table S1, Supporting Information). All subjects consumed each of the four model foods in a randomized order during four sessions (random Latin-square design). These sessions were separated by at least 3 weeks of wash-out because of the model foods provided higher amounts of vitamin D 3 and lutein than usual foods. The study was approved by the regional committee on human experimentation (IDRCB 2017-A01996-47, France): Clinical trial registry NCT03413267 (ClinicalTrials.gov). Procedures followed were in accordance with the Declaration of Helsinki of 1975 as revised in 1983. Objectives and requirements of the study were fully explained to all participants before beginning the study, and written informed consent was obtained from each subject. Participant flow chart is presented in Figure 1.

Model Food Preparation

The model foods were formulated with wheat flour (Francine T45, Grands Moulins de Paris, Ivry Sur Seine, France), extruded pea flour (Sativa 32/100, Sotexpro, Bermericourt, France), powdered sugar (Saint-Louis Sucre, Paris, France), sunflower oil (Lesieur, Asnières-sur-Seine, France), standard pasteurized egg yolk and granulated pasteurized egg white powders (Liot, Pleumartin, France), and sterilized water. Powdered vitamin D 3 was supplied by Nutrilo GmbH (ref 05-052, Cuxhaven, Germany). Lutein was provided from Diet Horizon (ref 2779, Bordeaux, France) as capsules containing a suspension of about 20 mg lutein in sunflower oil; a homogeneous dispersion of lutein was produced from these capsules and its concentration was determined by spectrophotometry at 445 nm using lutein molecular extinction coefficient.

Four model foods (one liquid and three solids) enriched in both vitamin D 3 and lutein were produced according to the flowcharts described in Figure S1, Supporting Information. Model food composition is presented in Table 2. Kneading, whisking, mixing steps, and heat treatment for custard were operated with a kitchen robot Thermomix TM5 (Vorwerk, Wuppertal, Germany). Cooking of biscuit, sponge cake, and pudding was operated with a semi-professional convection oven (De Dietrich, Niederbronn-Les-Bains, France). On dry matter, the model food recipes were all composed of 17% proteins, 52% carbohydrates, and 30% lipids. Custard was prepared as a packed-powder that was rehydrated and heated before use. Pudding, sponge cake, and biscuit were stored at -20 °C until use and then defrosted at room temperature.

Confocal Microscopy and Image Analysis of the Model Foods

Model food slices (300 µm) were cut either on a cryotome at -20 °C for biscuit, or on a microtome for sponge cake and pudding, and transferred onto glass slides. Fluorescent probes (Fast Green and Red Nile) were applied to the slides and let to diffuse for 10 min. The two probes were mixed with custard before incubation for 10 min, and the mixture was applied directly on slides. Slides were covered with a spacer and a cover slide and were imaged using a LSM880 inverted confocal microscope (Carl Zeiss AG, Oberkochen, Germany). Images were analyzed using Image J software. The results obtained were then analyzed as described previously [21] by fitting the cumulative percentage of lipid droplets area to the Rosin-Rammler equation. [22] This enabled, using linear regression techniques, to extract two parameters to describe the distribution curves of the lipid droplet sizes, namely ×50, which is the median particle area (µm²) and b, a dimensionless constant, which represent the broadness of the distribution.

Postprandial Experiments

Because it has been shown that season is an important prediction factor for vitamin D 3 response variation, [23] the postprandial experiments were performed in winter, and more precisely between November 2017 and March 2018. Subjects were asked to eat a standardized dinner between 7 pm and 8 pm the day before the postprandial experiment and to abstain from any food or beverage consumption afterward-with the exception of water. After an overnight fast, subjects arrived at the Unité d'Exploration Nutritionnelle and a baseline blood sample was taken. Subjects then consumed the model food together with 200 mL of mineral water. Subjects were asked to consume the model food within 15 min (to diminish the variability due to different rates of intake and, thus, gastric emptying). No other food was permitted over the following 8 h. However, subjects were asked to drink 100 mL water every 2 h (up to 500 mL). Additional blood samples were taken at 2, 3, 4, 5, 6, and 8 h after meal consumption. Blood was collected in heparin-lithium tubes that were immediately placed on ice and covered with aluminum foil to avoid light exposure, then centrifuged at 1500 × g for 10 min at 4 °C.

Chylomicron Preparation

Chylomicrons were prepared from plasma samples as previously described. [24] Immediately after recovery, chylomicrons were stored at -80 °C until analysis.

Lipid Analysis

Cholesterol, phospholipids, and triglycerides (TG) in chylomicrons were analyzed as previously described. [25]

Model Food Lutein and Vitamin D 3 Extraction and Analysis

Model foods (≈100 mg) were homogenized in 800 µL of saturated NaCl solution with 1-mm-diameter ceramic balls in 2 mL Eppendorf tubes using a FastPrep-24 5G (MP Biomedicals, Illkirch, France).

Lutein was extracted following the Bligh and Dyer method [26] and using apo-8'-carotenal as an internal standard. The extract was evaporated to dry under argon and dissolved in 600 µL of methanol/dichloromethane (65/35, v/v). A final volume of 10 µL was used for HPLC analysis.

Vitamin D 3 was extracted by adding 1 volume of ethanol, which also contained ergocalciferol (vitamin D 2 ) as an internal standard, and vortexing for 30 s. After adding 2 volumes of hexane, the mixture was vortexed for 10 min and centrifuged at 500 × g for 10 min at 4 °C. The upper phase (containing the different vitamin D forms) was collected and the sample was extracted a second time with hexane following the same procedure. The hexane phases were pooled and completely evaporated under nitrogen gas. The dried residue was dissolved in 200 µL methanol/dichloromethane (65/35, v/v) and a final volume of 100 µL was used for HPLC analysis.

Reverse-phase HPLC was then performed as previously described for lutein [27] and vitamin D 3 [28] quantification.

Chylomicron Lutein and Vitamin D 3 Extraction and Analysis

Lutein and vitamin D 3 from chylomicron samples were extracted by the methods previously described. [27,28] The whole extracts were evaporated to dryness under argon and dissolved in 100 µL of methanol/dichloromethane (65/35, v/v). A final volume of 90 µL was used for reverse-phase HPLC analysis. Lutein was detected at 445 nm and was identified by retention time and spectral analysis (300-500 nm) compared with pure lutein standard. Vitamin D 3 was detected at 265 nm and was identified by retention time and spectral analysis (200-400 nm) compared with pure vitamin D 3 standard. The two compounds were quantified by comparing peak area with standard calibration curves and were corrected by extraction efficiency based on recovery of internal standard.

Calculation and Statistical Analysis

The trapezoidal approximation method was used to calculate area under the curve (AUC) of the postprandial plasma chylomicron TG, lutein, and vitamin D 3 concentration over 8 h, henceforth referred to as TG, lutein, or vitamin D 3 response. Lutein and vitamin D 3 responses were corrected for ingested lutein and vitamin D 3 . Uncorrected curves are presented in Figures S2 andS3, Supporting Information. Coefficient of variation (CV) of TG, lutein, and vitamin D 3 responses were calculated as the ratio of standard deviation compared to the mean. All tests were two-sided, with a type I error set at 0.05. Continuous data were expressed as mean ± standard error of mean. The assumption of normality was assessed by using the Shapiro-Wilk test. Random-effects models for repeated data were performed to compare the four ingested matrices considering the following fixed effects: i) matrix and for longitudinal analyses, time-point evaluation and time x matrix interaction, and ii) taking into account between and within participant variability (subject as random-effect). Furthermore, due to cross-over randomized design, period, matrix order, and carry-over were also considered as covariates. When appropriate, a Sidak's type I error correction was applied to perform multiple comparisons. The normality of residuals from these models was studied using the Shapiro-Wilk test. When appropriate, a logarithmic transformation was proposed to achieve the normality of dependent outcome.

Statistical dependence between two variables was assessed by Pearson's correlation coefficient (r). Statistical analyses were performed using Stata software (version13, StataCorp, College Station, USA) and R software. For all tests, a p-value <0.05 was considered significant.

Results

Model Food Characteristics

The differences in the production process and the water content resulted in four strongly different model foods in terms of macroscopic aspect (Figure 2A). The biscuit was a very dry and brittle product. The sponge cake was a soft and dry gel with a honeycomb structure, while the pudding was a wet and compact gel. Finally, the custard was a thick liquid.

Confocal microscopy images of the four model foods are presented in Figure 2A. A continuous protein network was observed in biscuit, sponge cake, and pudding, as suggested by the uniform green background. In regard to the custard, protein aggregation was suspected because of the small green particles of about 50 µm. With regards to the lipid phase, fat droplets had different shapes depending on the model foods: fat globules could be seen in the custard versus fat "flakes" (i.e., large spots) in the biscuit. Both droplets and flakes were characterized by their area. Lipid droplets appeared bigger in custard and biscuit compared to sponge cake and pudding, in which lipids were more diffused. The distribution curves of lipid droplets area confirmed the bigger size of lipid droplets in the biscuit (Figure 2B): this matrix presented the higher median droplet size, around 18 000 µm. However, the median droplet size measured for custard was surprisingly the lowest amongst the four model foods (median size around 30 µm). Despite large lipid droplets were observed on the confocal images, the custard was actually mainly composed of very small fat droplets that could not be observed on the image presented in Figure 2A because of their reduced size. Sponge cake and pudding presented intermediate median area of lipid droplets, around 3000 and 600 µm, respectively. Lastly, for the distribution spread measured by the b constant, custard was found to have the narrowest size distribution (as indicated by the highest b value, 1.08), followed by pudding (0.46), sponge cake (0.38), and then biscuit (0.25).

Chylomicron TG Response to the Four Different Model Foods

The mean TG concentration after consumption of the different model foods is shown in Figure 3A. All subjects had extremely low chylomicron TG concentrations (<0.109 ± 0.064 mmol L -1 ) at T0, which confirmed that they were at fast at the beginning of the experiment.

The maximal chylomicron TG concentrations after custard, pudding, sponge cake, and biscuit consumption were obtained 3 h after meal intake (Figure 3B). The chylomicron TG responses (i.e., 0-8 h AUC) ranged as follows: 3.63 ± 0.62, 2.77 ± 0.44, 2.07 ± 0.33, and 1.63 ± 0.25 mmol L -1 per hour for custard, sponge cake, pudding, and biscuit, respectively. Custard TG response was significantly higher than pudding and biscuit responses (+75.3% and +123%, respectively, p < 0.0001). Sponge cake TG response was significantly higher than biscuit response (+69.7%, p = 0.002) (Figure 3C). The CV of the TG response was 64.8% and the ratio between the highest and the lowest responder was ≈ 16.

Chylomicron Vitamin D 3 Response to the Four Different Model Foods

The mean vitamin D 3 concentration after consumption of the different model foods is shown in Figure 4A and in Figure S2, Supporting Information. Ten out of the 12 subjects had vitamin D 3 concentrations below our limit of detection in their fasting chylomicron fraction, and the remaining two subjects had extremely low vitamin D 3 fasting concentrations on one experimental day (0.265 ± 0.073 nmol L -1 ).

The maximal chylomicron vitamin D 3 concentrations were obtained 3 h after custard, pudding, sponge cake consumptions and 4 h after biscuit consumption (Figure 4B). The chylomicron vitamin D 3 responses ranged as follows: 24.29 ± 3.31, 23.19 ± 2.83, 21.98 ± 3.08, and 19.18 ± 1.69 nmol L -1 per hour for sponge cake, pudding, custard, and biscuit, respectively. Sponge cake vitamin D 3 response was significantly higher than biscuit response (+26.6%, p = 0.047; Figure 4C). The CV of the vitamin D 3 response was 43.2% and the ratio between the highest and the lowest responder was ≈7. Considering all model foods, the vitamin D 3 response correlated positively with chylomicron TG response (Pearson's r = 0.496, p = 0.0003). When model foods were considered separately, the relationship between vitamin D 3 and TG responses was still significant for sponge cake (Pearson's r = 0.740, p = 0.006; Table 3).

Chylomicron Lutein Response to the Four Different Model Foods

The mean lutein concentration after consumption of the different model foods is shown in Figure 5A and in Figure S3, Supporting Information. Subjects had low fasting lutein concentrations in the chylomicron fraction (<6.419 ± 3.974 nmol L -1 ).

The maximal chylomicron lutein concentrations after custard, pudding, sponge cake, and biscuit consumption were obtained 3 h, 3 h 30 min, 3 h, and 4 h after meal intake, respectively (Figure 5B). The chylomicron lutein responses ranged as follows: 125.63 ± 25.25, 108.90 ± 26.45, 96.63 ± 19.20, and 76.08 ± 15.72 nmol L -1 per hour for sponge cake, custard, biscuit, and pudding, respectively. Lutein response was not significantly different depending on the model food (Figure 5C). The CV of the lutein response was 74.9% and the ratio between the highest and the lowest responder was ≈69.

Considering all model foods, the lutein response correlated positively with the chylomicron TG response (Pearson's r = 0.324, p = 0.025) and the chylomicron vitamin D 3 response (Pearson's r = 0.619, p < 0.0001). When model foods were considered separately, the positive relationship between lutein and TG responses was still significant in the case of pudding (Pearson's r = 0.607, p = 0.038), and the relationship between vitamin D 3 and lutein responses was still significant for custard and pudding (Pearson's r = 0.877, p = 0.0002 and Pearson's r = 0.650, p = 0.022, respectively; Table 3).

Discussion

In order to assess the bioavailability of lutein and vitamin D 3 , we measured the postprandial chylomicron lutein and vitamin D 3 responses after consumption of four different model foods. This approach is commonly used for experimental assessment of both vitamin D 3 [28] and carotenoid [29,30] absorption and bioavailability, as chylomicrons mostly carry fat-soluble components from intestinal origin. Only men were enrolled in the study. One of the reasons for choosing groups of the same sex lies in fact that it may decrease the variability of responses between individuals. This is particularly true in postprandial studies providing a lipid challenge because lipid or lipid phytochemical responses are often lower in women than men due to a more efficient chylomicron-remnant clearance [31] in relation with visceral adipose tissue accumulation. [32] The food matrix has long been considered as a key factor modulating the absorption of fat-soluble micronutrients/phytomicroconstituents. Most of the studies performed so far simultaneously explored the effects of both food composition and food structure on fat-soluble micronutrient bioavailability. Indeed, in vitro, both carotenoid transfer from plant matrix to mixed micelles and carotenoid absorption significantly vary from one plant matrix to another. [33] This effect was also observed in vivo as broccoli and green peas were more effective sources of lutein than spinach in humans. [34] This food matrix effect is less clear for vitamin D 3 . [35] A food matrix containing fat such as cheese did not increase vitamin D 3 bioavailability compared to a supplement in water in 8 healthy adults. [36] Similarly, vitamin D 3 bioavailability remained unaffected by a mushroom matrix compared to a supplement in 30 healthy adults. [37] In contrast, specific food components such as fatty acids [38] or phytosterols [39] were shown to modulate lutein bioavailability in humans. Similarly vitamin D 3 absorption was modulated in vitro or in rodents by fatty acids and phytosterols, as well as proteins [40] (see, for review, ref. [15]).

Besides, the effect of food processing on fat-soluble micronutrients bioavailability was specifically investigated. Although food processing leads to a partial degradation of vitamins and phytochemicals due to thermal and oxidative conditions, it also largely promotes matrix disruption. This disruption can balance the negative effect of micronutrient degradation by increasing micronutrient intestinal absorption. This effect is particularly important for carotenoids. Data obtained in humans highlighted that the consumption of 300 g per day of chopped spinach during 4 days improved lutein response compared to the consumption of whole leaf spinach. [34] Similarly, the absorption efficiency of the xanthophyll 𝛽-cryptoxanthin was 1.8-fold higher from orange juice than from fresh oranges. [41] By contrast, another clinical trial in which subjects received, for 3 weeks, spinach products including whole leaf spinach with an almost intact food matrix, minced spinach with a partially disrupted matrix, and enzymatically liquefied spinach, evidenced no effect of food matrix on lutein response. [42] These discrepancies can be due to the fact that i) the sources of spinach used were not the same in the different studies, and ii) food processing can affect several aspects of the food matrix, including its composition and its structure, in a manner that highly depends on processing parameters.

To specifically study the effect of food matrix structure on micronutrient postprandial response, we designed four model foods corresponding to real foods (i.e., biscuit, sponge cake, pudding, and custard). These model foods were prepared with the same ingredients but displayed contrasted microstructures and macrostructures. The four different model food recipes only differed in their water content. As water was shown to have a minor effect on gastric emptying, [43] this should not have impacted nutrient digestion. Nutrients were assessed in the final products to check the effect of processing. In particular, moderate differences in terms of lipid contents were observed after baking with a maximal difference of 11% was observed between pudding and sponge cake, which contained 40.3 ± 0.7 and 35.8 ± 1.4 g of fat per portion, respectively. This is likely due to lipid leaching at high temperature. However, a previous study performed in healthy volunteers showed that no difference was observed for TG postprandial responses after the ingestion of meals containing 30, 40, or 50 g of fat. [44] Consequently, the differences observed in our study in terms of (micro)nutrient responses, which were corrected by the amounts of ingested micronutrients, are expected to be mainly due to model food structure.

The CV for TG (64.8%), vitamin D 3 (43.2%), and lutein (74.9%) responses are consistent with those obtained in previous studies. In a previous trial performed in 39 healthy men, the CV of responses in chylomicrons were 49% for TG, [45] 47% for vitamin D 3 , [28] and 75% for lutein. [29] respectively. This confirms the very high interindividual variability regarding the absorption efficiency of fat-soluble compounds.

Interestingly, when considering all model foods, vitamin D 3 response was positively and significantly correlated with chylomicron TG response (Table 3). This is consistent with previous published data. [28] Lutein response was also significantly and positively correlated with chylomicron TG responses, in agreement with previous results. [29] Finally, lutein response also correlated with chylomicron vitamin D 3 response. These correlations were also significant when considering data for some individual matrices. Considering the existence of low, normal, and high lipid [46] or fat-soluble micronutrient responders, [47] our study also evidenced that the high TG responders are also likely high vitamin D 3 /lutein responders. Our study further indicates that high responders for one fat-soluble micronutrient are likely to be high responders for another fat-soluble micronutrients.

Previous data showed that supramolecular organization and physical state of lipids in native and formulated food products can modulate both the digestibility and metabolism of dietary fat. In particular, emulsification in small lipid droplets leads to an enhancement of fat absorption. [START_REF] Michalski | Steering Committee of RMT LISTRAL[END_REF] This was confirmed in our study as the highest TG postprandial response was observed after custard consumption (smallest fat droplets) while the lowest TG postprandial response was observed after biscuit (biggest fat droplets). Indeed, TG response after biscuit consumption was 55% lower than after custard consumption. This decrease could also be linked to the biscuit matrix as its solid structure would slow the accessibility of the lipid droplets to enzymes. However, none of the model foods had an impact on the time of the maximal postprandial TG blood concentration, which was observed at 3 h. This indicates that the structure of the food matrix neither delayed the accessibility of TG to enzymes during the digestion process, nor fatty acid absorption and secretion into chylomicrons.

Interestingly, such reduction in absorption after the ingestion of the biscuit matrix was also observed for vitamin D 3 (-21% after biscuit compared to sponge cake). This decrease in vitamin D 3 absorption may be due to a vectorization of this vitamin in a lipid emulsion displaying large-size fat droplets, which is are less efficiently absorbed. However, the model food structure did not significantly modify the lutein response. This lack of effect can be due to the above-mentioned extremely high interindividual variability regarding lutein responses (ratio between the highest and the lowest responder =69 for lutein, while this ratio =7 for vitamin D 3 ). The difference observed between vitamin D 3 and lutein responses may also be related to the fact that the two components display very different chemical structures. Indeed, vitamin D 3 is a sterol-like molecule while lutein is an oxygenated tetraterpenoid. Such difference can impact the ability of the molecules to transfer to mixed micelles and to be absorbed and secreted by enterocytes.

Both vitamin D 3 and lutein maximal blood concentrations were delayed for about 1 h after biscuit intake compared to other matrices. This result is consistent with a recent study in minipigs aiming at comparing lutein response after pure lutein (i.e., without matrix) or spinach lutein ingestion. This study showed no difference in lutein responses, but the lutein peak was delayed from 2 to 6 h with spinach, [49] indicating that the food matrix can modulate micronutrient absorption rate.

Additionally, both vitamin D 3 and lutein maximal blood concentrations were also delayed for 1 h compared to TG after biscuit intake (4 h vs 3 h). Such differences in maximal postprandial blood concentrations of TG and of vitamin D 3 /lutein suggest that, compared to other matrices, biscuit induced either i) a delay in intestinal absorption of fat-soluble micronutrients compared to TG, ii) a different mechanism for the incorporation of TG and fat-soluble micronutrients into chylomicrons, or iii) a differential clearance for these compounds in chylomicron remnants. Further investigations are required to validate these assumptions.

Finally, it is worth mentioning that the decrease in TG and vitamin D 3 responses observed in this study with some matrices may be balanced by increased TG and vitamin D 3 in the chylomicron secretion following the next meal. Indeed, the intestine has the capacity to store some lipid molecules such as cholesterol and vitamin E from one meal and to secrete them during subsequent postprandial periods. [50] In conclusion, we showed that food structure can influence the absorption of fat-soluble (micro)nutrients in humans. In particular, it can lead to modification in amounts of absorbed TG and vitamin D 3 . Thus, food structure might be a technological lever to optimize the bioefficiency of vitamin D 3 -supplemented foods in a context of deficient populations. Moreover, this study highlighted that food structure can lead to a delay in both vitamin D 3 and lutein postprandial peaks, although the physiopathological consequences of such absorption rate modulation are not yet known. Overall, this study confirms the importance of further investigation to better understand the impact of matrix characteristics on food nutritional quality. This may be of particular interest to develop functional foods for population subgroups deficient or sub-deficient in fat-soluble vitamins, such as the elderly.

Figure 1 .

 1 Figure 1. Participant flow chart.

Figure 2 .

 2 Figure 2. Model food characteristics. A) Macroscopic and microscopic images of model foods. Confocal images are at magnification ×20 and lipids appear in red and proteins in green. B) Lipid size distribution of model foods. Vertical lines correspond to ×50, the median particle area of each distribution (µm²).

Figure 3 .

 3 Figure 3. Chylomicron triglyceride (TG) responses. A) Baseline-adjusted chylomicron TG concentration over 8 h after the consumption of model foods. Data are means ± SEM (n = 12). At each time, a different letter indicates a significant difference (p < 0.05). B) Time of maximum TG concentration after the consumption of model foods. Vertical bars represent medians; boxes represent 25th-75th percentiles and whiskers indicate the maximal and the minimal values. A different letter indicates a significant difference (p < 0.05). C) Postprandial chylomicron TG responses (0-8 h AUC) after the consumption of model foods. Data are means ± SEM (n = 12). A different letter indicates a significant difference (p < 0.05).

Figure 4 .

 4 Figure 4. Chylomicron vitamin D 3 responses. A) Baseline-adjusted chylomicron vitamin D 3 concentration over 8 h after the consumption of model foods corrected by vitamin D 3 intake. Data are means ± SEM (n = 12). At each time, a different letter indicates a significant difference (p < 0.05). B) Time of maximum vitamin D 3 concentration after the consumption of model foods. Vertical bars represent medians; boxes represent 25th-75th percentiles and whiskers indicate the maximal and the minimal values. A different letter indicates a significant difference (p < 0.05). C) Postprandial chylomicron vitamin D 3 responses (0-8 h AUC) after the consumption of model foods corrected by vitamin D 3 intake. Data are means ± SEM (n = 12). A different letter indicates a significant difference (p < 0.05).

Figure 5 .

 5 Figure 5. Chylomicron lutein responses. A) Baseline-adjusted chylomicron lutein concentration over 8 h after the consumption of model foods corrected by lutein intake. Data are means ± SEM (n = 12). At each time, a different letter indicates a significant difference (p < 0.05). B) Time of maximum lutein concentration after the consumption of model foods. Vertical bars represent medians; boxes represent 25th-75th percentiles and whiskers indicate the maximal and the minimal values. A different letter indicates a significant difference (p < 0.05). C) Postprandial chylomicron lutein responses (0-8 h AUC) after the consumption of model foods corrected by lutein intake. Data are means ± SEM (n = 12). A different letter indicates a significant difference (p < 0.05).

Table 1 .

 1 Baseline characteristics of healthy men.

	Parameter	Value
	Age [y]	22.6 ± 0.4
	Weight [kg]	71.2 ± 1.4
	BMI [kg m -2 ]	2 2 . 1 ± 0.4
	Glucose [mmol L -1 ]	4.22 ± 0.12
	Hemoglobin [g dL -1 ]	1 5 . 6 ± 0.3
	25-Hydroxyvitamin D [µg L -1 ]	22.7 ± 1.1

Table 2 . Model food recipes and nutritional composition.

 2

	Custard	Pudding	Sponge cake	Biscuit

2000228 (3 of 11)

Mol. Nutr. Food Res. 2020, 2000228 © 2020 Wiley-VCH GmbH

Mol. Nutr. Food Res. 2020, 2000228 © 2020 Wiley-VCH GmbH 2000228 (11 of 11)

Acknowledgements

B.G. and M.H. contributed equally to this work. This project was funded by a grant from the Institut Carnot Qualiment. The authors are grateful to Dr. Claire Dufour (UMR SQPOV, Avignon, France) for helpful suggestions on the manuscript.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Supporting Information

Supporting Information is available from the Wiley Online Library or from the author.

Author Contributions

Conflict of Interest

The authors declare no conflict of interest.


François Lavoué
email: francois.lavoue@univ-grenoble-alpes.fr

Olivier Coutant

Pierre Boué

Laura Pinzon-Rincon

Florent Brenguier

Romain Brossier

Philippe Dales

Meysam Rezaeifar

Christopher J Bean

Adriano Čarman

Jérôme Cavaliere

Claudio Chèze

John Chiarabba

Glenn Clin- Ton

Wayne C Cougoulat

Luigia Crawford

Tibor Cristiano

Ezio Cz- Ifra

Stefania D'alema

Romuald Danesi

Anke Daniel

Iva Dannowski

Anne Dasović

Jean-Xavier Deschamps

Cécile Dessa

Sven Doubre

Ethz-Sed Egdorf

Tomislav Electronics Lab

Kasper Fiket

Wolfgang Fischer

Florian Friederich

Sigward Fuchs

Domenico Funke

Aladino Giardini

Zoltán Govoni

Gidera Grá- Czer

Stefan Gröschl

Ben Heimers

Davorka Heit

Marijan Herak

Johann Herak

Dejan Huber

Petr Jarić

Yan Jedlička

Hélène Jia

Edi Jund

Stefan Kissling

Bernhard Klingen

Petr Klotz

Heidrun Kolínský

Michael Kopp

Josef Korn

Lothar Kotek

Krešo Kühne

Dietrich Kuk

Jürgen Lange

Sara Loos

Deny Lovati

Lucia Malengros

Christophe Margheriti

Xavier Maron

Marco Martin

Francesco Massa

Thomas Mazzarini

Laurent Meier

Irene Métral

Milena Molinari

Anna Moretti

Jurij Nardi

Anne Pahor

Catherine Paul

Daniel Péquegnat

Damiano Petersen

Davide Pesaresi

Claudia Pic- Cinini

Thomas Piromallo

Jaroslava Plenefisch

Sil- Via Plomerová

Snježan Pondrelli

Roman Prevolnik

Marc Racine

Miriam Régnier

Joachim Reiss

Georg Ritter

Simone Rümpker

Marco Salimbeni

Werner Santulin

Sven Scherer

Detlef Schippkus

Vesna Schulte-Kortnack

Stefano Šipka

Daniele Solarino

Kathrin Spallarossa

Josip Spieker

Angelo Stipčević

Bálint Strollo

Gyöngyvér Süle

Eszter Szanyi

Christine Szűcs

Martin Thomas

Frederik Thorwart

Stefan Tilmann

Massimiliano Ueding

Luděk Valloc- Chia

René Vecsey

Joachim Voigt

Zoltán Wassermann

Chris- Tian Wéber

Viktor Weidle

Gauthier Wesztergom

Stefan Weyland

Felix Wiemer

David Wolf

Thomas Wolyniec

Mladen Zieke

Helena Žlebčíková Živčić

Understanding seismic waves generated by train traffic via modelling: implications for seismic imaging and monitoring

Trains are now recognized as powerful sources for seismic interferometry based on noise correlation, but the optimal use of these signals still requires a better understanding of their source mechanisms. Here we present a simple approach for modelling train-generated signals inspired by early work in the engineering community, assuming that seismic waves are emitted by sleepers regularly spaced along the railway and excited by passing train wheels. Our modelling reproduces well seismological observations of tremor-like emergent signals and of their harmonic spectra. We illustrate how these spectra are modulated by wheel spacing, and how their high-frequency content is controlled by the distribution of axle loads over the rail, which mainly depends on ground stiffness beneath the railway. This is summarized as a simple rule of thumb that predicts the frequency bands where most of train-radiated energy is expected, as a function of train speed and of axle distance within bogies. Furthermore, we identify two end-member mechanisms -single stationary source vs. single moving load-that explain two types of documented observations, characterized by different spectral signatures related to train speed and either wagon length or sleeper spacing. In view of using train-generated signals for seismic applications, an important conclusion is that the frequency content of the signals is dominated by highfrequency harmonics and not by fundamental modes of vibrations. Consequently, most train traffic worldwide is expected to generate signals with a significant high-frequency content, in particular in the case of trains travelling at variable speeds which produce truly broadband signals. By proposing a framework for predicting train-generated seismic wavefields over meters to kilometers distance from railways, this work paves the way for high-resolution passive seismic imaging and monitoring at different scales with applications to near-surface surveys (aquifers, civil engineering), natural resources exploration, and natural hazard studies (landslides, earthquakes and volcanoes).

Figures

Introduction

Trains have recently been recognized as powerful seismic noise sources for imaging (e.g. [START_REF] Nakata | Shear wave imaging from traffic noise using seismic interferometry by cross-coherence[END_REF][START_REF] Quiros | Seismic interferometry of railroad induced ground motions: body and surface wave imaging[END_REF][START_REF] Dales | Virtual sources of body-waves from noise correlations in a mineral exploration context[END_REF] and monitoring [START_REF] Brenguier | Train traffic as a powerful noise source for monitoring active faults with seismic interferometry[END_REF] the Earth's subsurface, from tens of meters to kilometers depth. By correlating train-generated noise between sensors, [START_REF] Nakata | Shear wave imaging from traffic noise using seismic interferometry by cross-coherence[END_REF], [START_REF] Quiros | Seismic interferometry of railroad induced ground motions: body and surface wave imaging[END_REF] and [START_REF] Dales | Virtual sources of body-waves from noise correlations in a mineral exploration context[END_REF] were able to recover body waves that can be processed to reveal subsurface structures, in a similar fashion as from active seismic data. [START_REF] Brenguier | Train traffic as a powerful noise source for monitoring active faults with seismic interferometry[END_REF] showed that train-generated body waves can be used to illuminate fault zones at a few kilometers depth, thereby allowing in situ monitoring of real faults in the field in similar conditions as those achieved in laboratory experiments (e.g. [START_REF] Scuderi | Precursory changes in seismic velocity for the spectrum of earthquake failure modes[END_REF]. However, the appropriate use of train-generated seismic noise for extracting information about the subsurface still requires a better understanding of their source mechanisms and of their spectral characteristics in the frequency band [1 -50] Hz of interest for seismological applications targeting the shallow crust.

While there exists detailed modelling studies in the engineering community of how trains generate ground vibrations in the vicinity of the railway (up to 100 m), including all the complexity of the interactions between trains, rail tracks and ballast (e.g. [START_REF] Kouroussis | Free field vibrations caused by highspeed lines: Measurement and time domain simulation[END_REF][START_REF] Li | Finite element model of ballasted railway with infinite boundaries considering effects of moving train loads and Rayleigh waves[END_REF], see also [START_REF] Connolly | Benchmarking railway vibrations -track, vehicle, ground and building effects[END_REF] for a review), there is a relative lack of modelling for understanding distant seismic observations [START_REF] Fuchs | Equidistant spectral lines in train vibrations[END_REF]Li et al., 2018a;[START_REF] Inbal | Sources of long-range anthropogenic noise in Southern California and implications for tectonic tremor detection[END_REF][START_REF] Brenguier | Train traffic as a powerful noise source for monitoring active faults with seismic interferometry[END_REF]. [START_REF] Fuchs | Equidistant spectral lines in train vibrations[END_REF] provided many clues for interpreting the signals recorded in a range of distances from the railway (300 to 1500 m). A striking observation is the presence of equidistant spectral lines in the spectrograms, which they explained as being caused by the repeated excitation of stationary sources (which they did not precisely identify) by the successive wheels of the train. They proposed that the fundamental frequency of the resulting harmonic signal may be expressed as f 1 = V train /∆ b and controlled by the train speed V train and the bogie distance ∆ b . Their explanations, however, are mainly based on observations and on rules of thumb derived from short-distance modelling [START_REF] Connolly | Benchmarking railway vibrations -track, vehicle, ground and building effects[END_REF][START_REF] Degrande | Free field vibrations during the passage of a Thalys high-speed train at variable speed[END_REF][START_REF] Ditzel | Seismograms of moving trains: comparison of theory and measurements[END_REF][START_REF] Wu | Vibration analysis of railway track with multiple wheels on the rail[END_REF][START_REF] Sheng | A comparison of a theoretical model for quasi-statically and dynamically induced environmental vibration from trains with measurements[END_REF]. Moreover, while these propositions seem to well explain most of the observations, they are not entirely satisfactory for explaining the larger spacing between frequency peaks or the Doppler effects that are sometimes observed (e.g. Fuchs et al., 2018, figs 6a,b and S2b). In this study, we further investigate these hypotheses via the simulation of seismic wave propagation.

The interest of this paper is threefold. First, we present an effective strategy for modelling train-generated signals. Without entering into the complexity of train-track interactions, our simple modelling strategy well reproduces the main characteristics of the signals recorded by seismic stations located at moderate distances from the railway (a few hundreds of meters to a few kilometers). Second, we clarify the mechanisms at play in the generation of seismic waves by trains. Although they are well known and documented in the engineering community (e.g. [START_REF] Krylov | Calculation of low-frequency ground vibrations from railway trains[END_REF], we feel a need to illustrate these mechanisms in relation to seismological observations. Third, we identify the main parameters that control these mechanisms, and therefore the amplitude and frequency content of the observed signals, and we highlight perspectives for seismic applications.

The outline of the paper is the following: we start by presenting our modelling strategy, based on 3D visco-elastic simulations and on Krylov's model of the quasi-static excitation generated by trains [START_REF] Krylov | Calculation of low-frequency ground vibrations from railway trains[END_REF]. We detail the design of the source time functions representing the train load applied on the ground and how we expect these functions and their spectral characteristics to control the final results. We then focus on two end-member cases: single moving load vs. single stationary source, and show that they generate spectral signatures that are very similar to the two categories of observations made by [START_REF] Fuchs | Equidistant spectral lines in train vibrations[END_REF]. Then we show how to reproduce signals similar to the observations, both in the time and frequency domain, when considering realistic cases of entire trains exciting many sources along the railway. We conclude by discussing the parameters that control the amplitudes and frequency content of train-generated signals, and we derive guidelines for detecting these signals in seismic ambient noise.

Modelling strategy

General approach and workflow

Following [START_REF] Krylov | Calculation of low-frequency ground vibrations from railway trains[END_REF] and [START_REF] Li | Finite element model of ballasted railway with infinite boundaries considering effects of moving train loads and Rayleigh waves[END_REF], we assume that the train load is applied onto the ground surface as point-source, stationary, quasi-static vertical forces acting on the sleepers, located every ∆ sleeper along the rail way, and activated by the passage of the train wheels (Fig. 1). Before going further, we shall specify that we use the term stationary in the same sense of spatially stationary as [START_REF] Fuchs | Equidistant spectral lines in train vibrations[END_REF], meaning that the sleepers are located at fixed positions and therefore constitute stationary sources. Besides, we describe the excitation as quasi-static because we only consider the effect of quasi-static pressure of wheel axles onto the track and ground [START_REF] Krylov | Calculation of low-frequency ground vibrations from railway trains[END_REF], and disregard dynamic excitations between wheels and track due for instance to wheel or rail irregularities (e.g. [START_REF] Sheng | A comparison of a theoretical model for quasi-statically and dynamically induced environmental vibration from trains with measurements[END_REF][START_REF] Lombaert | Ground-borne vibration due to static and dynamic axle loads of InterCity and high-speed trains[END_REF]Galvín et al., 2010a,b;[START_REF] Xia | Theoretical modeling and characteristic analysis of moving-train induced ground vibrations[END_REF]. Finally, the sleepers are small compared to the propagating wavelengths and are assumed to be rigid. Therefore we consider only one series of sources, indicated as a single rail in Fig. 1(b), which amounts to assume that the load of each train axle is applied as a single point source (i.e. the two wheels of each axle act simultaneously on each sleeper). Using source-receiver reciprocity, we simulate the signal that would be generated by a train passing over the sleepers and recorded by a sensor away from the railway as the sum of the signals resulting from a virtual source located at the sensor position and recorded at the sleepers positions (Fig. 1b).

Our workflow is therefore the following:

1. Simulation of wave propagation between one (virtual) source and all (virtual) receivers, using a Dirac delta source time function such as to convolve the resulting impulse responses with any desired source time function at the post-processing stage.

2. Lowpass filter (using a zerophase Butterworth filter) below the maximum frequency allowed by the simulation (f simulation max = 100 Hz).

3. Resample at 250 Hz (> 2f simulation max in order to satisfy Nyquist criterion). This is just to facilitate post-processing, given that the original simulation time step is usually of the order of 0.2 ms (5000 Hz) to satisfy the Courant-Friedrichs-Lewy (CFL) condition and ensure the stability of the numerical scheme.

4. Take time derivative to convert the simulated displacement (in m) into ground velocity (in m/s), which is the quantity that is usually recorded by seismic stations.

Note that this time derivative has an effect on the high-frequency content of the signals.

5. (optional) Select or mute specific arrivals (e.g. direct P or S, refracted P or S, surface waves).

6. Convolve individual impulse responses with source time functions representing the train passage over each sleeper.

7. Sum all the individual convolved seismograms to get the final seismogram resulting from the contributions of all sleepers.

Steps 6 and 7 can be formalized as

v z (x r , t) = N S j=1 S j (t) * G zz (x r , x j , t), (1)
where v z (x r , t) is the vertical component of the ground velocity recorded at receiver location x r along time t, N S is the number of sources (sleepers) along the railway, and * denotes time convolution between the source time function S j (t) of source j and the Green function G zz (x r , x j , t) of the medium between the location x j of the vertical source j and the receiver location x r .

In the following, we give more details on some stages of this workflow, in particular on the computation of the individual impulse responses (Green functions G zz ) and on the design of the source time functions S j (t) representing the train passage over each sleeper.

Computation of Green functions

The method used to compute the individual impulse responses before convolution and stacking is not critical to our argument and could be performed via a number of different techniques, either (semi)analytical or numerical. In this work we compute these individual impulse responses using the SEM46 software [START_REF] Trinh | Efficient timedomain 3D elastic and viscoelastic full-waveform inversion using a spectral-element method on flexible Cartesian-based mesh[END_REF] based on time-domain spectral elements under the visco-elastic approximation and for a maximum frequency of 100 Hz, so as to investigate the generation of high frequencies.

The use of spectral-element simulations presents numerous advantages. For instance, spectral elements enable an accurate localisation of the very densely-spaced virtual receivers and a consistent extraction of the simulated wavefield at these locations, on the same basis of Gauss-Lebatto-Legendre (GLL) polynoms as used in the modelling scheme. Second, the use of source-receiver reciprocity with a very high number ( 20,000) of virtual receivers (actual sources) partly mitigates the computational cost of this fullynumerical method. Needless to say, we only need a short simulation duration (4 s, i.e. the time for waves to travel from the source to the furthest receiver), and we do not simulate the entire duration of the train passage ( 540 s in the examples shown hereafter). Finally, spectral elements enable us to consider media with realistic physical properties (including attenuation) and with arbitrary heterogeneity. In this paper, we restrict ourselves to examples in a homogeneous medium, because we want to focus on the source and it turns out that the main mechanisms can be easily exemplified in such simple models (see Table 1 for the properties considered in our numerical experiments). Further work may investigate the effects of medium heterogeneity in real-scale cases, e.g. in southern California [START_REF] Brenguier | Train traffic as a powerful noise source for monitoring active faults with seismic interferometry[END_REF] or in our test sites in Ontario (Canada, [START_REF] Dales | Virtual sources of body-waves from noise correlations in a mineral exploration context[END_REF] and in the Dublin basin (Ireland, Rezaeifar et al., 2020).

Train source time functions

Characteristic lengths and expected frequencies at play

According to the considered mechanism consisting of sleepers hit by wheels, we may expect two fundamental frequencies, related 1. to the excitation of each source (sleeper) by the successive wheel pairs, spaced every ∆, resulting in a frequency f 1 = V train /∆, where the spacing ∆ may a priori be • the distance ∆ a between axles within a given bogie,

• the distance ∆ b 1 between two bogies of the same wagon (as suggested by [START_REF] Fuchs | Equidistant spectral lines in train vibrations[END_REF],

• the distance ∆ b 2 between two bogies of consecutive wagons (which is typically shorter than ∆ b 1 ),

• or the wagon length L w (as indicated by [START_REF] Krylov | Calculation of low-frequency ground vibrations from railway trains[END_REF].

2. to the passage of each wheel pair over successive sleepers, resulting in a frequency f 2 = V train /∆ sleeper called passage frequency by [START_REF] Krylov | Calculation of low-frequency ground vibrations from railway trains[END_REF]. Note that in this case, we consider the wheels as moving sources which may induce a Doppler effect.

Figure 1 depicts these characteristic lengths within a schematic train and Table 1 indicates the values considered in our numerical experiments. As a first task of this paper, let us define the source time function describing the passage of a train at a given source (sleeper) location, and clarify which spacing ∆ controls the frequency f 1 .

Effect of train geometry on frequency content

To define the source time function for a given sleeper, we must consider the passage of each wheel of the train over this particular sleeper. Let us start by considering a simplistic train with only one wheel axle per wagon, i.e. where the only characteristic distance is the wagon length, assumed the same for all wagons. In this case, a simple source time function can be defined as a Dirac comb with regular peaks at the arrival times of each wheel at the source location, with a period L w /V train related to wagon length L w and train speed V train . Figure 2 shows such a simple Dirac comb (in greeny yellow), together with its spectrum, corresponding to a 8-wagon train with 1 wheel pair per wagon hitting a single sleeper. The spectrum is harmonic, presenting peaks at frequencies that are multiples of the fundamental frequency f 1 = V train /L w . Because the train has a finite length (a finite number of wagons), the spectrum also presents some frequency leakage which causes the frequency peaks to broaden and not to be as sharp as their time-domain counterparts (this may have its importance later on, when convolving this Dirac comb with more sophisticated functions that might amplify minor peaks). We shall now ask ourselves what happens if we simulate a more realistic train, with 4 wheel axles per wagon (Fig. 1a).

Figure 2 presents in blue the Dirac comb that corresponds to a 8-wagon train with 4 axles per wagon, following the geometry of Fig. 1(a) and the values of Table 1. Intuitively, the time intervals between peaks reflect the distances between the wheels of the train. Less intuitively, we observe in the spectrum that adding more wheels (i.e. more characteristic lengths, see Fig. 1) does not introduce extra fundamental frequencies but modulates the original spectrum, as expected from the Fourier transform of a doublet of Dirac peaks (e.g. Krylov and Ferguson, 1994, eq. 18). This is shown in Fig. 2 as black dots for the modulation due to the axle distance ∆ a , whose maxima correspond to frequency f a = V train /∆ a (pink dashed line), and as black dashes for the modulation due to both the axle distance ∆ a and the bogie distance ∆ b 1 within one wagon, with peaks corresponding to frequency f b 1 = V train /∆ b 1 (purple dashed line). Although this description does not fully correspond to the train geometry (as it does not account for the bogie distance ∆ b 2 between consecutive wagons), note how the black dots and dashes well matches the amplitude of the blue peaks, while cancelling out some of the initial greeny-yellow peaks. As a side note, [START_REF] Krylov | Calculation of low-frequency ground vibrations from railway trains[END_REF] suggested that the zeroes introduced by the modulation effect may serve to suppress vibrations at chosen frequencies, for example at the passage frequency, by designing train and track geometry accordingly.

Another conclusion is that it is the wagon length L w that controls the main fundamental frequency f 1 , and not the interval between bogies as intuited by [START_REF] Fuchs | Equidistant spectral lines in train vibrations[END_REF]. This leads us to suggest that the high-speed trains generating the signals observed by [START_REF] Fuchs | Equidistant spectral lines in train vibrations[END_REF] at AlpArray station A002A, with a spectral line spacing of ∆f = f 1 = 1.25 Hz, are not travelling at a speed of 85 km/h, as the authors suggested based on bogie distance, but rather at 119.25 km/h ( 120 km/h), which is therefore the train speed we will consider in our numerical experiments.

These simple Dirac combs provide us with a first understanding of the source time functions and of their spectral characteristics. However, they are not entirely satisfactory, as they assume that the load of each wheel is applied as an instantaneous point source on a single sleeper. As a consequence, assuming simple Dirac combs tends to over-estimate the high-frequency content of the source time functions, and therefore of the resulting seismograms. For more realistic source time functions, we need to account for the fact that the load of each axle is actually distributed over several sleepers, via the elastic response of the rail track (e.g. [START_REF] Krylov | Calculation of low-frequency ground vibrations from railway trains[END_REF][START_REF] Li | Finite element model of ballasted railway with infinite boundaries considering effects of moving train loads and Rayleigh waves[END_REF].

Effect of the distribution of axle loads over sleepers on frequency content

Figure 3 presents examples of functions which may be used to describe the distribution of the load of a single wheel over rail track and sleepers. The first approximation discussed in the previous section is a simple Dirac comb (in blue). More realistically, we may consider that, due to the elasticity of the track, point-source loads like wheels actually affect several sleepers simultaneously. A first way to take this into account is to consider simple weights as given for instance by Paderno (2009, their fig. 4) on 5 sleepers (black line in Fig. 3). Considering a load distributed over a certain width has an effect on the spectrum of the temporal loading function and tends to decrease its high-frequency content (intuitively, the wider the distribution, the lower the frequency content). We can expect to retrieve this effect in the final signal that will result from the convolution of the Green function with a source time function (such as the ones shown in Fig. 2) that will itself be convolved with this loading function.

Alternatively, we may consider a smoothed version of these weights, for instance using a Gaussian function (in green in Fig. 3). It is important here to notice that not only the width but also the shape of the loading function is critical in controlling the highfrequency content of the spectrum: for a comparable width, a Gaussian function will yield much less high frequencies (> 20 Hz) than a triangular function.

The use of ad hoc loading functions is therefore dubious and may result in incorrect results in the high-frequency band that is of interest to us. As a consequence, we base our source time functions on a more physical approach and consider as a loading function the (opposite of the) counteracting force exerted by the track on the wheel, which is also the force exerted by the track on the ground via the sleepers, as given by the Euler-Bernoulli elastic beam (E-BEB) model (Krylov and Ferguson, 1994, eqs 2 and 7;Li et al., 2018a, eqs 2 and 3):

P j (x i ) = 2F i e -β|x i -x j | cos(β|x i -x j |) + sin(β|x i -x j |) ∆ sleeper x 0 . (2)
In this formulation, the loading function P j (x i ) corresponds to the vertical force (in N) at sleeper j caused by the axle load of wheel i, and depends on the distance x i -x j between wheel i and sleeper j, on the axle load F i (assumed known for a given train weight), on the sleeper spacing ∆ sleeper (quite standard and assumed to be ∆ sleeper = 24" = 0.6096 m in our numerical experiments, see Table 1), on the elastic modulus E (in N/m 2 ) and cross-sectional momentum I (in m 4 ) of the rail (assumed known for standard steel rails), and finally on the stiffness α (in N/m 2 ) of the ground beneath the rail, via the coefficient

β = (α/4EI) 0.25 (in m -1
) and the total deflection distance x 0 = π/β (in m, [START_REF] Krylov | Calculation of low-frequency ground vibrations from railway trains[END_REF]. Among these parameters, we identify the ground stiffness α as the main parameter susceptible to vary in real-life situations, depending not only on the actual stiffness of the ground locally, but more generally on the coupling between the rail track, the sleepers and the ground, via the ballast and substratum. Figure 3 shows loading functions given by the E-BEB model with various values of ground stiffness: [START_REF] Li | Finite element model of ballasted railway with infinite boundaries considering effects of moving train loads and Rayleigh waves[END_REF], in greeny yellow), α = 61.8 MN/m 2 (in orange, after Krylov and Ferguson, 1994, who document this as a typical value for British railways), and α = 800 MN/m 2 (in red, which matches the weights of Paderno, 2009 and -more importantly-best matches the seismic observations, as we will see later in Fig. 4b,c). The stiffer the ground beneath the rail is, the narrower is the spatial loading function, and the higher the frequencies in the source time functions, and therefore -we expect-in the final seismograms. Moreover, amplitudes are normalised in Fig. 3. In absolute, ground stiffness also has an effect on the magnitude of the force exerted on the ground: the stiffer the ground beneath the rail is, the larger the vertical force exerted by the track on the ground (see Fig. S1).

α = 10 MN/m 2 (
In conclusion, our final source time function S j (t) (in orange in Fig. 2) is the result of the convolution of the Dirac comb corresponding to 4 wheel axles per wagon (in blue in Fig. 2) with the time derivative of the loading function P j (t) given by the E-BEB model:

S j (t) = dP j (t) dt * N axles i=1 δ(t -t ij ), (3)
where t ij is the arrival time of wheel i at sleeper j, and the time derivative is needed to simulate ground velocities instead of displacements. In summary, just looking in detail at the design and characteristics of this source time function already enables us to draw the following conclusions (or, more precisely, to retrieve the conclusions of [START_REF] Krylov | Calculation of low-frequency ground vibrations from railway trains[END_REF]:

1. The sharp frequency peaks observed by [START_REF] Fuchs | Equidistant spectral lines in train vibrations[END_REF] correspond to a fundamental frequency f 1 = V train /L w controlled by train speed and wagon length L w , which is the true spatial periodicity of the train wheels.

2. The spectrum of the source time functions is modulated by the spacing of the wheels within each wagon, in particular by the distance ∆ a between axles within each bogie.

3. The frequency bandwidth of the source time functions is further controlled by the distribution of the load of each wheel axle over the rail track and sleepers, which can be estimated through the E-BEB model and parameterised through an effective ground stiffness (which in practice may also reflect the coupling between track and ground): the stiffer the ground (the better the coupling), the narrower the loading function, and the higher the frequency content.

4. Whether we look at ground displacement or velocity matters since the time derivative relating both also has an effect on the spectra.

5. From points 2, 3 and 4, we may derive a rough rule of thumb according to which we expect most of the energy of ground velocity spectra to be in the frequency band [ 0.5 f a , 1.5 f a ], where f a = V train /∆ a is controlled by the train speed and the distance ∆ a between axles within each bogie. Secondary energy packets can also be expected in subsequent frequency bands of the form [ (k -0.5) f a , (k + 0.5) f a ], in between zero-energy notches located at frequencies (k + 0.5) f a . This is illustrated in Fig. 2 S3).

We expect these characteristics of the source time functions to be reflected in the final seismograms generated by train passages, which we shall now investigate in more details by adding the contributions of all sleepers.

Modelling results

Two end-member cases: single stationary source vs. single moving load 2018) interpreted the frequency spacing as ∆ f 2 = 2 f 1 , i.e. as being related to f 1 = V train /L w and therefore controlled by wagon length L w (Fuchs et al., 2018, figs 6a,b), while the Doppler frequency shift rather suggests the effect of a moving load. Yet, the frequency spacing ∆ f 2 = 2.54 Hz suggests that the characteristic length at play is not as short as the sleeper spacing, but may be a longer spatial period related to irregularities of the railway (e.g. unwelded rails). For instance, considering that loads are applied only every 19 sleeper introduces a spatial period of 19∆ sleeper 11.6 m, therefore resulting in a fundamental frequency f 2 = f 2 /19 = 2.52 Hz (Fig. 4f) which would reasonably explain the observations (Fig. 4e).

These two end-member cases, single stationary source vs. single moving load, therefore enable us to understand the two main mechanisms at play, and nicely relate to the two types of observations made by [START_REF] Fuchs | Equidistant spectral lines in train vibrations[END_REF], at least in terms of their spectral characteristics. But they are only end-member, simplistic cases, and the question now is: what happens when considering all wheels exciting all sleepers? How do we generate signals similar to the observations, both in the frequency domain (harmonic spectra) and in the time domain (tremor-like emergent shapes)?

Realistic cases (all sleepers, all wheels): effect of sleeper regularity Figure 4(f) shows the signal obtained with a realistic wheel geometry (see Fig. 1a) and loads applied every 19 sleeper along the railway, with a regular spacing of 19∆ sleeper 11.6 m, for a train with 6x 23-m-long wagons moving at 105 km/h (such as to try and match the observations in Fig. 4(e), which [START_REF] Fuchs | Equidistant spectral lines in train vibrations[END_REF] attributed to a local commuter train, see their table 1)1 . Surprisingly, the spectrogram is much closer to the case of a single moving load (Fig. 4d, see also Fig. S4c-d), with a clear Doppler effect and a spectrum dominated by frequency f 2 = V train /(19∆ sleeper ), than to the stationary source case (Fig. 4a). However, this result matches well the few observations made by [START_REF] Fuchs | Equidistant spectral lines in train vibrations[END_REF] (their figs 6a, 6b and S2b) where the Doppler effect is visible. In this specific case, it is indeed easy to get confused and to interpret the main spectral line spacing in terms of multiples of f 1 but a close inspection of the spectrum indicates that the frequency peaks are much better explained by the moving load mechanism and the frequency f 2 /19 than by frequency f 1 , which barely plays a role here (further tests in this sense are presented in the Supplemental Material, Fig. S4).

The question is then: how can we reproduce seismograms similar to Fig. 4(b), with a tremor-like emergent signal in the time-domain and a spectrum dominated by frequency f 1 = V train /L w , controlled by wagon length L w ?

Figure 4(c) shows the signal obtained with all wheels and all sleepers, with slightly irregular sleepers spaced every ∆ sleeper = 0.6096 ± 0.05 m (i.e. we add a random perturbation to the original spacing, with a maximum perturbation of ±5 cm, which [very roughly] corresponds to a sleeper half-width). Again, we consider a train with 8x 26.5-mlong wagons, moving at 119.25 120 km/h such that f 1 = 1.25 Hz, as observed by [START_REF] Fuchs | Equidistant spectral lines in train vibrations[END_REF]. A striking observation here is that it is sufficient to assume slightly irregular sleepers to reproduce a signal and spectrum that are very close to most of the observations made by [START_REF] Fuchs | Equidistant spectral lines in train vibrations[END_REF] on station A002A (emergent tremor-like time-domain signal and harmonic spectrum dominated by frequency f 1 = V train /L w , controlled by wagon length L w ). While it does not dominate the spectrum anymore, the peak at frequency f 2 = V train /∆ sleeper still appears because of the mean sleeper spacing ∆ sleeper = 0.6096 m. This f 2 -peak is not visible in Fig. 4(b) because it is above the Nyquist frequency of 50 Hz but we expect that it should be seen in seismograms resulting from trains travelling at slower speeds. We actually suspect that it is responsible for the secondary maximum that is sometimes observed around 40 Hz by [START_REF] Fuchs | Equidistant spectral lines in train vibrations[END_REF] in the spectra generated by commuter trains (e.g. Fuchs et al., 2018, fig. 3). Finally, we speculate that our modelling parameter of sleeper regularity that controls the predominance of one or the other endmember mechanisms (single stationary source vs. single moving load) actually reflects the properties of the railway infrastructure in real cases, i.e. not only the actual regularity of sleeper locations but also their coupling with the ground via the ballast and substratum.

Discussion

Limitations of the modelling strategy

A careful reader will have noticed that our modelling does not quantitatively reproduce the absolute amplitudes of observed ground velocities: while the ground velocities generated by high-speed trains observed at station A002A are of the order of 10 -6 to 10 -5 m/s (Fig. 4b,e), the signals simulated in similar configurations have amplitudes of the order of 10 -8 m/s (Fig. 4c,f). This discrepancy may be partly related to poorly-controlled parameters such as sensor-ground coupling or possible site effects, or to our choice of medium properties, which consist of low quality factors and high seismic velocities in order to reduce computation costs, while in reality station A002A is located in a sedimentary basin where soft, low-velocity materials may amplify the amplitude of surface waves. Future work will focus on quantitative comparisons with data acquired in better-controlled environments in order to calibrate our synthetic amplitudes.

It may also be argued that this discrepancy between observed and simulated amplitudes could be due to the fact that our modelling does not include all the mechanisms involved in reality. Indeed, by following Krylov's model [START_REF] Krylov | Calculation of low-frequency ground vibrations from railway trains[END_REF], we only consider a quasi-static excitation caused by the axle load on the rail track, and we neglect dynamic effects related to other excitation mechanisms due for instance to wheel and rail roughness, rail joints, and sleeper-ground coupling (e.g. [START_REF] Xia | Theoretical modeling and characteristic analysis of moving-train induced ground vibrations[END_REF]. Similarly, we do not account for the generation of horizontal forces due to friction and rail curvature, for the heavier axle loads of locomotives at the front and end of the train, or for the elasticity and spatial extension of the sleepers which actually support two rails while we considered only one (which, again, amounts to assume that the sleepers are rigid and that the load of each train axle is applied as a single point source). In hindsight, we think that the discrepancy between observed and simulated amplitudes is most likely due to the latter point, i.e. to the fact that we inject energy via point sources in the spectral-element simulations, without considering the surface of the sleepers.

Nevertheless, we have shown that our simple modelling considering only quasi-static excitations is able to reproduce -qualitatively at least -the seismic observations, which is also consistent with [START_REF] Krylov | Calculation of low-frequency ground vibrations from railway trains[END_REF]'s statement that the quasi-static approximation is sufficient in the frequency range [0 -100] Hz. We are confident that this modelling strategy is sufficient to catch the main mechanisms at play and enables us to investigate the general behaviour of generated amplitudes with respect to some parameters. Because we do not reproduce absolute amplitudes accurately, however, we restrict the discussion of amplitudes to a qualitative analysis in the following, and present related results in Supplemental Material (Figs S2 andS3).

Parameters controlling signal amplitudes: implications for seismic detection

In view of using train-generated signals for seismic applications, it is important to assess in which conditions we can expect to detect these signals. Having all the above assumptions and limitations in mind, we may therefore attempt to understand the parameters that control the amplitudes of train-generated signals. According to our approach and assumptions, signal amplitude depends on

• the stiffness of the ground beneath the rail: the magnitude of the vertical force acting on the sleepers increase non-linearly with ground stiffness (except for its quasi-static component which depends only on axle load, see Fig. S1).

• the length and total weight of the train, with two different effects:

via the weight of individual wagons which will reflect on the axle load F i carried by each wheel: eq. 2 shows that the magnitude of the vertical force acting on the sleepers increase linearly with axle load F i , and so will the amplitude of the resulting signal.

via the number of wagons, as the amplitude of the signals also increase with the number of wheel loads N wl . More precisely, our numerical experiments show that the maximum amplitude of the time-domain signals increases as √ N wl , while the energy in a given frequency bandwidth increases linearly with the number of wheel axles (see Fig. S2).

• the train speed: faster trains make more noise. More precisely, our numerical experiments show that both the maximum time-domain amplitude and the total signal energy increase supra-linearly with train speed, and so does the energy in frequency bands that evolve proportionally with train speeds, such as [ 0.5 , 1.5 ] f a , with f a = V train /∆ a (see Fig. S3). Interestingly, frequency bands of the form [ k -0.5 , k + 0.5 ] f a always contain approximately the same proportion of the spectrum, which brings us back to our rule of thumb stating that most of the energy of train-generated signals is expected to be in the [ 0.5 , 1.5 ] f a frequency band controlled by train speed and axle distance ∆ a within each bogie. A direct implication of this for detecting train-generated signals in seismic records is that the considered frequency band should be adapted to the estimated speed of the trains generating these signals.

In conclusion, if we want to maximize the detection of train-generated signals, we should targets portions of railways that enable a good transmission of vibrations to the ground (i.e. with sleepers well coupled to a stiff ground) and where long and/or heavy trains travel at relatively high, and ideally variable, speeds. Moreover, the considered frequency band should be adapted to the estimated speed of the trains generating the signals, targeting the frequency band [ 0.5 , 1.5 ] f a where we expect most of the energy, or higher-frequency bands of the form [ k -0.5 , k + 0.5 ] f a , and avoiding narrow frequency bands around the zeroes of the modulated spectrum ([k + 0.5] f a and, in a lesser extent, [k + 0.5] f b 1 ). If the trains travel at variable speeds, their signal becomes more broadband and the effect of these notches might be less critical (see next section). Of course, the detection of these signals in field recordings will also depend on the surrounding noise level, as well as on medium heterogeneity, attenuation and scattering.

Variable train speed

Figure 5 presents the example of a signal generated by a train travelling at a variable speed. When the speed of the train varies, the fundamental frequencies f 1 = V train /L w and f 2 = V train /∆ sleeper vary too, as well as their harmonics, which results in undulating spectral lines in the spectrograms, and in a blurring of the frequency peaks (and notches) in the corresponding spectra, where it is then difficult to identify individual rays. The spectrum is not harmonic anymore but becomes truly broadband as frequency notches disappear, which is actually an advantage for seismic applications based on noise crosscorrelations, as distinct spectral rays are difficult to handle by seismic interferometry. Note however that our rule of thumb remains useful: the frequency band [ 0.5 , 1.5 ] f a , where f a can be estimated from a maximum or averaged train speed, still contains most of the energy of the signal. Note also that the variable speed might lead to other prominent maxima in the spectrum, as for instance around 37 Hz in Fig. 5 where a change of train speed happens to occur when the train is closest to the station, which enhances the secondary maximum related to the Doppler effect.

Wave propagation effects

So far we have eluded wave propagation effects on purpose in order to focus on source mechanisms. Of course, once generated, train vibrations propagate as any other seismic wave. For the sake of completeness, we discuss here how wave propagation effects impact train-generated signals. We present associated results in Supplemental Material (Figs S5 andS6).

Short-vs. long-distance signals

Our tests suggest that long-range signals recorded at several kilometers distance have the same properties as the examples we have shown here at short range (300 to 1200 m). In particular, the main characteristics of the spectra obtained at short vs. long distances are nearly identical. The only differences between short-and long-distance signals that we could notice are related to foreseeable wave propagation effects, namely (1) to the expected decay of high frequencies at long distances, and (2) to the relative decay of amplitudes: short-distance signals have a shorter apparent duration because their relative amplitudes are dominated by the effect of the few closest sleepers, while long-distance signals have a longer apparent duration because they see a larger number of sources contributing with an equivalent importance (Figs S5a,b).

Effect of seismic properties (velocities and attenuation)

Obviously, intrinsic attenuation also has an effect on the relative decay of amplitudes and of high frequencies in the spectra. Besides, higher velocities also result in a longer apparent duration of the signal, which can be explained in terms of relative decay of amplitudes due to geometrical spreading over a smaller number of longer wavelengths. Less intuitively, higher velocities yield lower amplitudes, due to the fact that the amplitude of surface waves is enhanced in soft, low-velocity media. As already mentioned, this may be one of the reasons why the amplitudes of observed train-generated signals are larger than our synthetic amplitudes, obtained for seismic velocities that have been chosen artificially high to mitigate computation costs. Finally, a last effect of seismic velocities relates to the Doppler frequency shift, which will increase if seismic velocities get closer to the (apparent) train speed (i.e. if seismic velocities decrease, or if train speed increases, or if the sensor is aligned with the train trajectory).

Wave partitioning (body vs. surface waves)

We also conducted tests in a medium with a velocity gradient where we selected either the refracted P wave (first arrival) or the surface waves (most energetic arrivals) in the impulse response, before convolving with the train source time functions and stacking. These tests aimed at mimicking long-distance wave propagation in heterogeneous media where highfrequency surface waves would vanish due to strong scattering and attenuation in the shallow subsurface, which is difficult to reproduce in our synthetic simulations. Again, we observe nearly identical spectra in both cases, which suggests that our conclusions, although mostly based on short-distance signals dominated by surface waves, remain valid for long-distance body waves that sample deeper regions of the crust [START_REF] Brenguier | Train traffic as a powerful noise source for monitoring active faults with seismic interferometry[END_REF].

Finally, it is worth mentioning that simulations performed in similar configurations but in media with homogeneous vs. gradient velocity structures yield very similar results, suggesting that medium structure plays little role in the patterns of the resulting signals, which complexity mostly comes from the source, i.e. from the train source time functions and from the interferences between waves generated by the many sources involved.

Conclusions

In this paper, we have presented a modelling strategy for simulating seismic waves generated by trains, explained the main mechanisms at play, and derived implications for seismic applications. Our results showed that simple assumptions following Krylov's model of a quasi-static excitation due to axle loads applied as vertical forces distributed on discrete sleepers are sufficient to generate signals that reproduce very well the main characteristics of seismic observations (emergent tremor-like signals with harmonic spectra). Our modelling enabled a detailed understanding of the spectral characteristics of train-generated signals and of their fundamental frequencies. As already noted by [START_REF] Krylov | Calculation of low-frequency ground vibrations from railway trains[END_REF], the exact frequency content of the signals results from a modulation of the initial harmonic spectrum by the wheel spacing within each wagon, and their high-frequency content is further controlled by the distribution of the load of each wheel over the track and sleepers, which most notably depends on the stiffness of the ground beneath the railway (i.e. on track-ground coupling).

Moreover, we have identified two end-member mechanisms that correspond to the two types of observations made by [START_REF] Fuchs | Equidistant spectral lines in train vibrations[END_REF] and depend on sleeper regularity (i.e. on the properties of the railway infrastructure in real cases). In the case of perfectly regular sleepers, signals are dominated by the signature of each wheel acting as a single moving load passing over the sleepers. The spectrum then displays peaks corresponding to the fundamental frequency f 2 = V train /∆ sleeper controlled by train speed and sleeper spacing ∆ sleeper , and to its harmonics. In the case of irregular sleepers, signals are dominated by the response of each sleeper acting as a stationary source excited by the passage of the wheels. The spectrum then displays peaks at harmonics of the fundamental frequency f 1 = V train /L w controlled by train speed and wagon length L w .

Beyond these technical explanations, we were able to derive important consequences for the future use of train-generated signals for seismic applications. First, it is worth noting that trains generate a very broad and high frequency content [1 -50 Hz or above], as observed by (e.g.) [START_REF] Fuchs | Equidistant spectral lines in train vibrations[END_REF][START_REF] Brenguier | Train traffic as a powerful noise source for monitoring active faults with seismic interferometry[END_REF]; [START_REF] Dales | Virtual sources of body-waves from noise correlations in a mineral exploration context[END_REF], because most of the energy actually comes from high-frequency harmonics and not from the fundamental frequencies, which are not necessarily very high (f 1 0.25 to 5 Hz, depending on train speed and wagon length). To summarize this, we have derived a simple rule of thumb according to which most of the energy radiated by trains should be expected in the range [ 0.5 , 1.5 ] f a , where the frequency f a = V train /∆ a is controlled by train speed and by the distance ∆ a between axles within a bogie (usually of the order of ∆ a 1.5 to 3 m). Higher frequency bands of the form [ k -0.5 , k + 1.5 ] f a are also expected to contain some energy which may prove particularly useful when looking specifically for body waves (e.g. [START_REF] Nakata | Body wave extraction and tomography at Long Beach, California, with ambient-noise interferometry[END_REF][START_REF] Brenguier | Train traffic as a powerful noise source for monitoring active faults with seismic interferometry[END_REF][START_REF] Dales | Virtual sources of body-waves from noise correlations in a mineral exploration context[END_REF]. Second, we were able to identify the parameters that control the amplitudes of train-generated signals, and therefore to give guidelines in order to maximize their detection: an ideal seismic source is a long, heavy and fast train traveling on a railway whose sleepers are well coupled to a stiff soil.

Of course, the detection of train-generated signals will also depend on the noise level due to natural or other anthropic noise sources, as well as on medium heterogeneity, attenuation and scattering. While we cannot give general and quantitative rules, we believe that this paper provides elements of understanding and guidelines that will help future case-specific studies to detect and process train-generated seismic signals in ambient noise data (e.g. [START_REF] Nakata | Shear wave imaging from traffic noise using seismic interferometry by cross-coherence[END_REF][START_REF] Quiros | Seismic interferometry of railroad induced ground motions: body and surface wave imaging[END_REF][START_REF] Inbal | Sources of long-range anthropogenic noise in Southern California and implications for tectonic tremor detection[END_REF]Li et al., 2018a;[START_REF] Dales | Virtual sources of body-waves from noise correlations in a mineral exploration context[END_REF], or to design experiments that will explicitly target these signals (e.g. [START_REF] Nakata | Shear wave imaging from traffic noise using seismic interferometry by cross-coherence[END_REF][START_REF] Quiros | Seismic interferometry of railroad induced ground motions: body and surface wave imaging[END_REF][START_REF] Brenguier | Train traffic as a powerful noise source for monitoring active faults with seismic interferometry[END_REF].

Data and resources

Simulations have been performed with the SEM46 code developed in the SEISCOPE Consortium (http://seiscope2.osug.fr). Most of the computations were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is partly supported by the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche. This study is based on data from the AlpArray Seismic Network (2015) which at the time of publication was not publicly available. Please visit the project page (http://www. alparray.ethz.ch/en/seismic_network/backbone/data-policy-and-citation) for more details on data access. All data processing and plotting were done using Python and the ObsPy toolbox [START_REF] Krischer | ObsPy: A bridge for seismology into the scientific Python ecosystem[END_REF].

Supplemental Material for this article includes (1) more details on the effect of ground stiffness on the amplitude of the loading force according to the E-BEB model, (2) examples of signals dominated by the signature of moving loads, (3) details on the effect of train length and train speed on the amplitude of the signals discussed in the paper, and ( 4) examples of wave propagation effects mentioned in the discussion. For the sake of reproducibility, the computer programs developed for the convolutions of Green functions with train source time functions are available at https://gricad-gitlab.univ-grenoblealpes.fr/pacific/publications/2020_Lavoue-et-al_SRL_supplemental-material.

All websites were last accessed on September 2020.

Tables

Medium properties V P = 3400 m/s, V S = 2000 m/s, ρ = 2600 kg/m 3 , Q P = 100, Q S = 50 ground stiffness beneath railway: α = 800 MN.m -2 1). Note that seismic velocities are deliberately quite high in order to reduce computational costs. Attenuation is also deliberately high in order to better mimic a finite-size medium. We have checked that our main conclusions do not depend on medium properties.

Train properties

L w = 26.5 m, ∆ a = 2.5 m, ∆ b 1 = 19 m, ∆ b 2 = 7.5 m or L w = 23 m, ∆ a = 2.3 m, ∆ b 1 = 18.5 m, ∆ b 2 = 4.

  (b) where most of the energy of the final source time function (in orange) lies between 7 and 20 Hz, while another packet of energy lies between 20 Hz and 33 Hz. It is further illustrated in the Supplemental Material(Figs S2 and

Figure 4

 4 Figure 4(a) shows the seismogram (top), spectrogram (bottom left) and spectrum (right) obtained with a single stationary source (one sleeper only) excited by all wheels of a train with 8x 26.5-m-long wagons, moving at 119.25 km/h ( 120 km/h), as suggested by revisiting Fuchs et al. (2018)'s calculation (see above). Note the zoom on the time axis which makes visible the individual spikes corresponding to the signal of each wheel pair. Since only one sleeper is involved, the spectrum is very similar to the spectrum of the individual source time functions (Fig.2), with a spectral line spacing related to frequency f 1 = V train /L w , controlled by the wagon length L w , and a modulation related to frequency f a = V train /∆ a , linked to axle distance ∆ a . We indicate frequency f 2 =

  Figures

Figure 1 :

 1 Figure 1: Geometry of a train (a) and of the numerical experiments (b).

Figure 2 :

 2 Figure2: Source time functions (a) and corresponding spectra (b) for a sleeper located at x = 1 km along the railway, consisting of pure Dirac combs for 1 and 4 wheel axles per wagon (in greeny-yellow and blue, respectively) and of a 4-axle Dirac comb convolved with a loading function derived from the E-BEB model, assuming a ground stiffness α = 800 MN/m 2 (in red, see also Fig.3). The orange curves correspond to the final source time function obtained after time differentiation of the red curve, in order to simulate ground velocities instead of displacements. We therefore expect the orange spectrum in (b) to be representative of the energy content of the final seismogram obtained by summing the contributions of all sources. In panel (b), black dots and dashes represent the modulation due to axle distance and bogie distance within wagons, respectively. Associated fundamental frequencies are indicated as vertical dashed lines.

Figure 3 :

 3 Figure 3: Loading functions representing the distribution of a single wheel load on the track, as a function of distance between the wheel and a given sleeper (a) and corresponding spectra (b), after conversion of the space axis into time using a train speed of 120 km/h.

Figure 4 :

 4 Figure 4: Seismograms (top of each panel), spectrograms (bottom of each panel) and spectra (right of each panel) resulting from train passages. The first column (panels a,d) corresponds to the two synthetic end-members, the second column (panels b,e) to observations, and the third column to realistic synthetic cases involving all sleepers and all train wheels. The first row (panels a, b, c) corresponds to signals dominated by the end-member signature of a single stationary source, while the second row (panels d,e,f) presents signals dominated by the signature of a single moving load. AlpArray data (b,e) are highpass-filtered above 0.5 Hz. Synthetic data are lowpass-filtered below 100 Hz. Spectrograms are computed over time windows of 5 s for observations (b,e) and 20 s for synthetics (a,c,d,f) with a 90% overlap. Synthetic spectra are computed on tapered data (blue seismograms in panels c, d, f) in order to avoid the start and the end of the signal where only some of the train wheels are involved (in orange).

Figure 5 :

 5 Figure 5: Seismogram (top), spectrogram(bottom) and spectrum (right) resulting from a train travelling over slightly irregular sleepers at a speed increasing from 50 km/h to 120 km/h between 0 and 120 s, then decreasing to 95 km/h ( 230 s), increasing again to 120 km/h ( 320 s) and finally decreasing to 50 km/h. The spectrogram is computed over time windows of 5 s with a 90% overlap. The spectrum is computed on tapered data (blue seismograms) in order to avoid the start and the end of the signal where only some of the train wheels are involved (in orange in the upper panel).

Table 1 :

 1 Acquisition and railway geometry sleeper spacing ∆ sleeper = 24" = 0.6096 m number of sleepers N sleeper = 19358 railway length L railway = 11800 m 12 km distance between sensor and railway d min = 300 to 1200 m Settings of the numerical experiments presented in the paper. Train properties are taken fromFuchs et al. (2018, Table

	5 m

Here our choice of train parameters aims at presenting the example that best illustrates the underlying mechanisms but we realize that it makes the comparison between Figs4(a,c) and Figs4(d,f) not straightforward. For the sake of completeness, we present more results in the Supplemental Material (Fig.S4) for one-to-one comparison.

Acknowledgments

This research emanates from PACIFIC -Passive seismic techniques for environmentally friendly and cost-effective mineral exploration -which has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No. 776622 (https://www.pacific-h2020.eu). We also acknowledge support from the European Research Council under grant No. 817803, FAULTSCAN. ISTerre is part of Labex OSUG@2020 (ANR10 LABX56).

We thank the Editors, an anonymous Reviewer, and Florian Fuchs for their review of the manuscript. We further thank Florian Fuchs for providing us with the timing information of the seismograms used in his paper (Fuchs et al., 2018) and for fruitful discussions at the occasion of the Cargèse workshop on Passive Imaging 2019. We acknowledge the sponsors of the SEISCOPE Consortium in which the SEM46 code has been developed. We thank Ludovic Métivier for his advice on the software and Marwan Irnaka for the use of the Seiscope PyTools.

We acknowledge the use of data from the AlpArray network and would like to thank the AlpArray Seismic Network Team: György HETÉNYI, Rafael ABREU, Ivo ALLE-GRETTI, Maria-Theresia APOLONER, Coralie AUBERT, Simon BESANÇON, Maxime BÈS DE BERC, Götz BOKELMANN, Didier BRUNEL, Marco CAPELLO, Martina