Olivier Mazet

Willy Rodríguez

Lounès Chikhi
email: lounes.chikhi@univ-tlse3.fr

Mazet Olivier

Rodríguez Valcarce

Chikhi Lounès

Demographic inference using genetic data from a single individual: separating population size variation from population structure

Keywords: symmetric island model, population size change, maximum likelihood estimation, demographic history, coalescence time

INTRODUCTION

The sheer amount of genomic data that is becoming available for many organisms with the rapid development of sequencing technologies represents new opportunities for population genetics research. It is hoped that genomic data will increase our ability to reconstruct the history of populations [START_REF] Li | Inference of human population history from individual wholegenome sequences[END_REF] and detect, identify and quantify selection (Vitti et al. 2013). While this increase in genetic information will likely help biologists and anthropologists to reconstruct the demographic history of populations, it also exposes old challenges in the field of population genetics. In particular, it becomes increasingly necessary to understand how genetic data observed in present-day populations are influenced by a variety of factors such as population size changes, population structure and gene flow [START_REF] Nielsen | Statistical inferences in phylogeography[END_REF]. Indeed, the use of genomic data does not necessary lead to an improvement of statistical inference. If the model assumed to make statistical inference is fundamentally mis-specified, then increasing the amount of data will lead to increased precision for perhaps misleading if not meaningless parameters and will not reveal new insights [START_REF] Nielsen | Statistical inferences in phylogeography[END_REF][START_REF] Chikhi | The confounding effects of population structure, genetic diversity and the sampling scheme on the detection and quantification of population size changes[END_REF][START_REF] Heller | The confounding effect of population structure on Bayesian skyline plot inferences of demographic history[END_REF].

For instance, several recent studies have shown that the genealogy of genes sampled from a deme in an island model is similar to that of genes sampled from a non structured isolated population submitted to a demographic bottleneck [START_REF] Peter | Distinguishing between population bottleneck and population subdivision by a Bayesian model choice procedure[END_REF][START_REF] Chikhi | The confounding effects of population structure, genetic diversity and the sampling scheme on the detection and quantification of population size changes[END_REF][START_REF] Heller | The confounding effect of population structure on Bayesian skyline plot inferences of demographic history[END_REF]. As a consequence, using a model of population size change for a spatially structured population may falsely lead to the inference of major population size changes (Städler et al. 2009;[START_REF] Peter | Distinguishing between population bottleneck and population subdivision by a Bayesian model choice procedure[END_REF][START_REF] Chikhi | The confounding effects of population structure, genetic diversity and the sampling scheme on the detection and quantification of population size changes[END_REF][START_REF] Heller | The confounding effect of population structure on Bayesian skyline plot inferences of demographic history[END_REF][START_REF] Paz-Vinas | The demographic history of populations experiencing asymmetric gene flow: combining simulated and empirical data[END_REF]. Conversely, assuming a structured model to estimate rates of gene flow when a population has been submitted to a population size change, may also generate misleading conclusions, even though the latter case has been much less documented. More generally, previous studies have shown that spatial processes can mimic selection [START_REF] Currat | Comment on "Ongoing adaptive evolution of ASPM, a brain size determinant in Homo sapiens" and "Microcephalin, a gene regulating brain size, continues to evolve adaptively in humans[END_REF], population size changes [START_REF] Chikhi | The confounding effects of population structure, genetic diversity and the sampling scheme on the detection and quantification of population size changes[END_REF][START_REF] Heller | The confounding effect of population structure on Bayesian skyline plot inferences of demographic history[END_REF] or that changes in gene flow patterns can mimic changes in population size (Wakeley 1999;[START_REF] Broquet | Genetic bottlenecks driven by population disconnection[END_REF]. The fact that such dissimilar processes can generate similar coalescent trees poses exciting challenges [START_REF] Nielsen | Statistical inferences in phylogeography[END_REF]. One key issue here is that it may be crucial to identify the kind of model (or family of models) that should be used before estimating and interpreting parameters.

One solution to this problem is to identify the "best" model among a set of competing models. This research program has been facilitated by the development of approximate Bayesian computation (ABC) methods [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF][START_REF] Beaumont | Approximate Bayesian computation in population genetics[END_REF][START_REF] Cornuet | Inferring population history with DIY ABC: a user-friendly approach to approximate Bayesian computation[END_REF][START_REF] Beaumont | Approximate Bayesian computation in evolution and ecology[END_REF]. For instance, using an ABC approach, [START_REF] Peter | Distinguishing between population bottleneck and population subdivision by a Bayesian model choice procedure[END_REF] showed that data sets produced under population structure can be discriminated from those produced under a population size change by using up to two hundred microsatellite loci genotyped for 25 individual. In some cases, relatively few loci may be sufficient to identify the most likely model [START_REF] Sousa | May)Population divergence with or without admixture: selecting models using an ABC approach[END_REF][START_REF] Peter | Distinguishing between population bottleneck and population subdivision by a Bayesian model choice procedure[END_REF]), but in others, tens or hundreds of loci may be necessary [START_REF] Peter | Distinguishing between population bottleneck and population subdivision by a Bayesian model choice procedure[END_REF]. ABC approaches are thus potentially very powerful but it may still be important to improve our understanding of the coalescent under structured models.

In the present study we are interested in describing the properties of the coalescent under two models of population size change and population structure, respectively, and in devising a new statistical test and estimation procedures. More specifically we re-derive the full distribution of T 2 , the time to the most recent common ancestor for a sample of size two for a model of sudden population size change and for the n-island model. We then use a maximum likelihood approach to estimate the parameters of interest for each model (timing and ratio of population size change former and rate of gene flow and deme size for the latter). We develop a statistical test that identifies data sets generated under the two models. Finally, we discuss how these results may apply to genomic data and how they could be extended to real data sets (since T 2 is not usually known) and other demographic models.

In particular we discuss how our results are relevant in the context of the PSMC (Pairwise Sequentially Markovian Coalescent) method [START_REF] Li | Inference of human population history from individual wholegenome sequences[END_REF], which has been now

Demographic models

Population size change: We consider a simple model of population size change, where N (t) represents the population size (N , in units of genes or haploid genomes) as a function of time (t) expressed in generations scaled by N , the population size, and where t = 0 is the present, and positive values represent the past (Figure 1 (a)). More specifically we assume a sudden change in population size at time T in the past, where N changes instantaneously by a factor α. This can be summarized as

N (t) = N (0) = N 0 for t ∈ [0, T [, N (t) = N (T ) = αN 0 for t ∈ [T, +∞[. If α > 1 the
population went through a bottleneck (Figure 1) whereas if α < 1 it expanded. Since N represents the population size in terms of haploid genomes, the number of individuals will therefore be N/2 for diploid species. Note also that for a population of constant size the expected coalescence time of two genes is N generations, which therefore corresponds to t = 1. We call this model the SSPSC, which stands for Single

Step Population Size Change.

[Figure 1

The distribution of coalescence times: qualitative and quantitative analyses

In this section we used previous results [START_REF] Herbots | Stochastic models in population genetics: genealogy and genetic differentiation in structured populations[END_REF][START_REF] Donnelly | Coalescents and genealogical structure under neutrality[END_REF] to derive the distribution of coalescent times for the two models of interest. We show that even though they are different, these distributions can be similar under an indefinitely large number of parameter values (Figures 2 and3). Moreover we show that even when the distributions are distinguishable, their first moments may not be. In particular, we show that the first two moments (mean and variance) are near identical for a large number of parameter combinations. Before doing that we start by providing a simple intuitive rationale explaining why and how a model of population structure can be mistaken for a model of population size change. This intuitive approach is important because it allows us to understand how the parameters of the two models ((T , α) and (M , n), respectively) are linked.

Intuitive and qualitative rationale: We start by taking two genes sampled in the present-day population under the Single Step Population Size Change (SSPSC) model. If we assume that α > 1 (population bottleneck from an ancient population of size N 1 to a current population of size N 0 , with N 1 = αN 0 ) the probability that the two genes coalesce will vary with time as a function of N 0 , N 1 and T . If T is very small, then most genes will coalesce at a rate determined by N 1 , whereas if T is very large the coalescence rate will be mostly determined by N 0 . If we now take two genes sampled from the same island in the Structured Symmetrical Island (StSI) model, we can also see that their coalescence rate will depend on N , the size of the island and on m, the migration rate. If m is very low, the coalescence rate should mostly depend on N . If m is high, the two genes may see their lineages in different islands before they coalesce. As a consequence the coalescence rate will depend on the whole set of islands and therefore on the product nN , where n is the total number of islands. This intuitive description suggests that there is an intrinsic relationship between T and 1/M , and between α and n. The reason why structured populations exhibit signals of bottlenecks is because in the recent past the coalescence rate depends on the local island size N , whereas in a more distant past it depends on nN . In other words, it is as if the population size had been reduced by a factor of n. As we will see this rationale is only qualitatively correct, but it suggests that if we want to distinguish them it may be necessary to derive the full distribution of the coalescence times under the two models. We shall denote these coalescence times T SSP SC 2 and T StSI 2 , respectively.

Derivation of the distribution of coalescence times:

The distribution of T SSP SC 2 : The generalisation of the coalescent in populations of variable size was first rigorously treated in Donnelly and Tavare (1995)), and is clearly exposed in Tavaré ( 2004)). If we denote by λ(t) the ratio N (t) N (0) where t is the time scaled by the number of genes (i.e. units of coalescence time, corresponding to N (0)t generations),

we can compute the probability density function (pdf ) f SSP SC T 2 (t) of the coalescence time

T SSP SC 2
of two genes sampled in the present-day population. Indeed, the probability that two genes will coalesce at a time greater than t is

P(T SSP SC 2 > t) = e -t 0 1 λ(x) dx , (1)
where

λ(x) = I [0,T [ (x) + αI [T,+∞[ (x),
and

I [a,b[ (x) is the Kronecker index such that I [a,b[ (x) =        1 for x ∈ [a, b[ 0 otherwise.
Given that the pdf is

f SSP SC T 2 (t) = (1 -P(T SSP SC 2 > t)) P(T SSP SC 2 > t) = e -t I [0,T [ + e -T -1 α (t-T ) I [T,+∞[ .
This leads to the following pdf

f SSP SC T 2 (t) = e -t I [0,T [ (t) + 1 α e -T -1 α (t-T ) I [T,+∞[ (t). ( 2
)
The distribution of T We can compute this probability as follows. Taking two genes from the same island and going back in time, there are three events that may occur: a coalescence event (with rate 1), a mutation event (with rate θ) and a migration event (with rate M ). Taking now two genes from different islands, they cannot coalesce and therefore only a mutation or a migration event may occur. Migration events can then bring the lineages in the same island with probability 1 n-1 , and in different islands with probability n-2 n-1 . We thus obtain the following coupled equations:

p s (θ) = 1 1 + M + θ + M 1 + M + θ p d (θ),
and

p d (θ) = M/(n -1) M + θ p s (θ) + M (n -2)/(n -1) M + θ p d (θ).
10

By solving them, we obtain

p s (θ) = θ + γ D and p d (θ) = γ D with γ = M n -1 and D = θ 2 + θ(1 + nγ) + γ.
We can then obtain the full distribution through the Laplace transform formula, if we note that

p s (θ) = θ + γ (θ + α)(θ + β) = a θ + α + 1 -a θ + β with a = γ -α β -α = 1 2 + 1 + (n -2)γ 2 √ ∆ , where α = 1 2 1 + nγ + √ ∆ and β = 1 2 1 + nγ - √ ∆ .
Noting now that for any θ and any α we have +∞ 0 e -αs e -θs ds = 1 θ + α , it is straightforward to see that the pdf of T StSI 2 is an exponential mixture: . Their analytic values can be easily expressed as functions of the model parameters:

f StSI T 2 (t) = ae -αt + (1 -a)e -βt . (3)
E T SSP SC 2 =1 + e -T (α -1), V ar T SSP SC 2 =1 + 2T e -T (α -1) + 2αe -T (α -1) -(α -1) 2 e -2T , E(T StSI 2 ) =n, V ar(T StSI 2 ) =n 2 + 2(n -1) 2 M .
It is interesting to note that the expected time in the StSI model is n and does not depend on the migration rate [START_REF] Durrett | Probability models for DNA sequence evolution[END_REF]. The variance is however, and as expected, a function of both n and M . For the SSPSC model, the expected coalescence time is a function of both T and α. We note that it is close to 1 when T is very large and to α when T is close to zero. Indeed, when the population size change is very ancient, even if α is very large the expected coalescence time will mostly depend on the present-day population size, N 0 . Similarly, when

T is small it will mostly depend on N 1 . The relationship that we mentioned above between n and α (and between M and 1/T ) can be seen by noting that when T is close to zero (and M is large), the expectations under the two models are α and n, and the variances are

V ar T SSP SC 2 ≈ 1 + 2α(α -1) -(α -1) 2 = α 2 and V ar(T StSI

2

) ≈ n 2 . This exemplifies the intuitive rationale presented above. This relationship is approximate and will be explored below, but can be illustrated in more general terms by identifying scenarios with similar moments.

As figure 2 shows, the two models provide near-identical pairs of values for (E(T 2 ), V ar(T 2 ))

for "well chosen" parameters (T, α) and (M, n). Here by setting T to 0.1 (and M to 9, i.e.

1/M ≈ 0.11) whereas α and n were allowed to vary from 1 to 100, and from 2 to 100, respectively, we see that the two models exhibit very similar behaviours. We also plotted a second example obtained by setting M to 0.5 and T to 1.09, and varying n and α as above.

These examples ilustrate how n and α (respectively, M and 1/T ) are intimately related.

[Figure 2 about here.]

The near-identical values obtained for the expectation and variance under the two models explains why it may be difficult to separate models of population size change from models of population structure when the number of independent genetic markers is limited. However, the differences between the distributions of coalescence times under the two models suggest that we can probably go further and identify one model from another. For instance, figure 3 shows that even in cases where the first two moments are near-identical (T = 0.1 and α = 10 versus M = 7 and n = 9), it should be theoretically possible to distinguish them.

This is exactly what we aim to do in the next section. In practice, we will assume that we have a sample of n L independent T 2 values (corresponding to n L independent loci ) and will use these T 2 values to (i) estimate the parameter values that best explain this empirical distribution under the two models of interest, (ii) use a statistical test to compare the empirical distribution with the expected distribution for the ML estimates and reject (or not) one or both of the models. For simplicity, and to make it easier to read, we will often use the term loci in the rest of the manuscript when we want to mention the number of independent T 2 values.

[Figure 3 about here.]

Model choice and parameter estimation

General principle and parameter combinations: Given a sample (t 1 , ..., t n L ) of n L independent observations of the random variable T 2 , we propose a parameter estimation procedure and a goodness-of-fit test to determine whether the observed distribution of the T 2 values is significantly different from that expected from the theoretical T SSP SC

2 or T StSI 2 distributions.
This sample can be seen as a set of T 2 values obtained or estimated from n L independent loci.

We took a Maximum Likelihood (ML) approach to estimate the parameters (T, α) and (M, n)

under the hypothesis that the n L -sample was generated under the T SSP SC ). Note that the estimation procedure and the KS test were performed on independent sets of T 2 values. We thus simulated twice as many T 2 values as needed (2n L instead of n L ). With real data that would require that half of the loci be used to estimate ( T , α) and ( M , n), whereas the other half would be used to perform the KS test.

We expect that if the estimation procedure is accurate and if the KS test is performing well we should reject the SSPSC (respectively, the StSI) model when the data were simulated under the StSI (resp., the SSPSC) model. On the contrary we should not reject data simulated under the SSPSC (resp., the StSI) model when they were indeed simulated under that model. To validate our approach we used (t 1 , ..., t 2n L ) data sampled from the two T 2 distributions and quantified how the estimation procedure and the KS test performed. In order to do that, we varied the parameter values ((T, α) and (M, n)) for various 2n L values as follows. For T and α we used all 36 pairwise combinations between these two sets of values (0.1, 0.2, 0.5, 1, 2, 5), and (2,4,10,20,50,100), respectively. For M and n we used all the 48 combinations between the following values (0.1, 0.2, 0.5, 1, 5, 10, 20, 50) and (2,4,10,20,50,100), respectively. For 2n L we used the following values (40,100,200,400,1000,2000,20000). Altogether we tested 588 combinations of parameters and number of loci. For each 2n L value and for each parameter combination (T, α) (or (M, n)) we realized 100 independent repetitions of the following process. We first simulated a sample of 2n L values using the pdf s of the SSPSC (resp. StSI) model with (T, α) (resp. (M, n)). We then used the first n L values to obtain the ML estimates ( T , α) for the SSPSC model and ( M , n)

for the StSI model. Then, we performed a KS test using a 0.05 threshold on the second half of the simulated data (i.e. n L values) with each of the theoretical distributions defined by the estimated parameters. Finally, after having repeated this process 100 times we recorded all estimated parameters and counted the number of times we rejected the SSPSC and StSI models for each parameter combination and each 2n L value.

Maximum likelihood estimation in the SSPSC case: We know from section the pdf of the coalescence time in the SSPSC model of two genes. We can thus write the likelihood function for any couple of parameters (α, T ), given one observation t i as:

L t i (α, T ) = 1 α e -T -1 α (t i -T ) I [0,t i [ (T ) + e -t i I ]t i ,+∞[ (T ).
Given n L independent values t = (t 1 , t 2 , ..., t n L ), the likelihood is:

L SSP SC (α, T ) = n L i=1 L t i (T, α),
and taking the log it gives:

log(L SSP SC (α, T )) = n L i=1 log(L t i (α, T )).
Lemma 0.1 Given a set of n L independent observations {t 1 , t 2 , ..., t n L }, the log-likelihood function log(L SSP SC ) has no critical points in R 2 .

For the proof and some comments, see Supplementary Materials.

As a consequence of this lemma, we take (α, T ) = argmax a∈{1,...,n L } {log(L SSP SC (m a ))} as the Maximum Likelihood estimates, where

m a =       n L i=1 t i I ta<t i -Kt a K + 1 , t a       , a ∈ {1, 2, ..., n L }. with K = n L i=1 I t i <ta
Maximum likelihood estimation in the StSI case: Under the StSI model the expression of the critical points is not analytically derived. We know from section the pdf of coalescence times for two genes. Given n L independent values t = (t 1 , t 2 , ..., t n L ) we can compute the log-likelihood function for any set of parameters (n, M ) as:

log(L StSI (n, M )) = n L i=1 log(ae -αt + (1 -a)e -βt ))
We used the Nelder-Mead method (Nelder and Mead 1965) implementation of scipy [START_REF] Jones | SciPy: Open source scientific tools for Python[END_REF]) to find numerically an approximation to the maximum of the likelihood function. This method returns a pair of real numbers (n, M ). Since n should be an integer we kept either n or n + 1, depending on which had the largest log-likelihood value. represents the amount of information (the number of T 2 values or independent loci.) The second thing to note is that the two parameters are very well estimated when we use 10, 000 values of T 2 . This is particularly obvious for n compared to α, probably because n must be an integer, whereas α is allowed to vary continuously. For instance, for most simulations we find the exact n value (without error) as soon as we have more than 1000 loci. However, we should be careful in drawing very general rules. Indeed, when fewer T 2 values are available (i.e. fewer independent loci), the estimation precision of both parameters depends also on T and M , respectively. Interestingly, the estimation of α and n are remarkable even when these parameters are small. This means that even "mild" bottlenecks may be very well quantified (see for instance the Supplementary materials for α = 2, T values between 0.1 and 1 when we use only 1000 loci). We should also note that when the bottleneck is very old (T = 5) the estimation of the parameters is rather poor and only starts to be reasonable and unbiased for n L = 10, 000. This is not surprising since the expected T M RCA is 1. Under the SSPSC model most genes will have coalesced by t = 5, and should therefore exhibit T 2 values sampled from a stationary population (i.e. α = 1). As the number of loci increases, a small proportion will not have coalesced yet and will then provide information on α. The expected proportion of genes that have coalesced by T = 5 is 0.993. the values shown here we can see that T , the age of the bottleneck is very well estimated even when α = 2 (for n L = 10, 000). In other words, even a limited bottleneck can be very precisely dated. For stronger bottlenecks fewer loci (between 500 and 1000) are needed to still reach a high precision. This is particularly striking given that studies suggest that it is hard to identify bottlenecks with low α values [START_REF] Girod | Inferring Population Decline and Expansion From Microsatellite Data: A Simulation-Based Evaluation of the Msvar Method[END_REF]. Interestingly, the panels (b), (d) and (f) seem to suggest that it may be more difficult to estimate M than T . As we noted above this observation should be taken with care. Indeed, T and M are not equivalent in the same way as α and n. This is why we chose to represent a value of M such that M = 1/T , and why one should be cautious in drawing general conclusions here. Altogether this and the previous figure show that it is possible to estimate with a high precision the parameters of the two models by using only 500 or 1000 loci from a single diploid individual. There are also parameter combinations for which much fewer loci could be sufficient (between 50 and 100).

In Figure 6 we show some results of the KS test for the two cases (See the Supplementary Materials for the other parameter combinations). In the left-hand panels ((a), (c), and (e))

the data were simulated under the SSPSC model and we used the StSI model as a reference (i.e. we ask whether we can reject the hypothesis that genetic data were generated under a structured model when they were actually generated under a model of population size change). In the right-hand panels ((b), (d) and (f)) the same data were compared using the SSPSC model as reference and we computed how often we rejected them using a 5% rejection threshold. The left-hand panels exhibit several important features. The first is that, with the exception of T 2 = 5 we were able to reject the wrong hypothesis in 100% of the cases when we used 10, 000 independent T 2 values.

This shows that our estimation procedure (as we saw above in figures 4 and 5) and the KS test are very powerful. The second feature is that for T = 5, the test performs badly whatever the number of independent loci (at least up to 10, 000). This is expected since the expected T M RCA of two genes is t = 1, and 99.3% of the loci will have coalesced by t = 5.

This means that out of the 10, 000, only c.a. 70 loci are actually informative regarding the pre-bottleneck population size. Another important feature of the left-hand panels is that the best results are generally obtained for T = 1, 0.5 and 2, whichever the value of α. This is in agreement with [START_REF] Girod | Inferring Population Decline and Expansion From Microsatellite Data: A Simulation-Based Evaluation of the Msvar Method[END_REF] in that very recent population size changes are difficult to detect and quantify. The observation is valid for ancient population size changes as well. The right-hand panels are nearly identical, whichever α value we used (see also Supplementary Materials), and whichever number of T 2 values we use. They all show that the KS test always rejects a rather constant proportion of data sets. This proportion varies between 3 and 15%, with a global average of 8.9%. Altogether our KS test seems to be conservative. This is expected because for low n L values the estimation of the parameters will tend to be poor. Since the KS test uses a reference distribution based on the estimated rather than the true values, it will reject the simulated data more often than the expected value of 5%. show, as in Figure 6, that for all the values of n L and n we reject a rather constant proportion of data sets (between 5 and 10%). Altogether the two previous figures (figures 6 and 7) show that it is possible to identify the model under which genetic data were generated by using genetic data from a single diploid individual. [Figure 4 about here.] effective population size changes do. Hence, we think that our work could be helpful to the aim of setting these two scenarios apart in order to detect (for example) false bottleneck signals. Nevertheless, while our study provides several new results, there are still several important issues that need to be discussed and much progress that can still be made.

T 2 and molecular data

The first thing to note is that we assume, throughout our study, that we have access to the coalescence times T 2 . In real data sets, this is never the case and the T 2 are rarely estimated from molecular data. While this may seem as a limitation, we should note that the recent method of [START_REF] Li | Inference of human population history from individual wholegenome sequences[END_REF], that uses the genome sequence of a single individual to infer the demographic history of the population it was sampled from, actually estimates the distribution of T 2 values. In its current implementation the PSMC software does not output this distribution but it could be modified to do it. Note however, that the PSMC should only be able to provides a discretized distribution in the form of a histogram with classes defined by the number of time periods for which population size estimates are computed.

In any case, this suggests that it is in theory possible to use the theoretical work of Li and

Durbin to generate T 2 distributions, which could then be used with our general approach.

Moreover, it should be possible to use the theory developed here to compute, conditional of the T 2 distribution, the distribution of several measures of molecular polymorphism. For instance, under an infinite site mutation model it is in principle possible to compute the distribution of the number of differences between pairs of non recombining sequences for the two demographic models analysed here. Similarly, assuming a single stepwise mutation model it should be possible to compute the distribution of the number of repeat differences between two alleles conditional on T 2 . However, we must add here that while it is easy to use these distributions to simulate genetic data (rather than T 2 values) it is not straightforward to then use the genetic data to estimate the models' parameters and apply a test to identify the model under which the data were simulated. This would probably require a Khi-2 test since discrete rather continuous distributions would be compared. This is an issue that would deserve a full and independent study.

Demographic models

In our study we limited ourselves to two simple but widely used models. It would thus be important to determine the extent to which our approach could apply to other demographic models. The n-island or StSI model is a widely used model and it was justified to use it here. One of its strongest assumptions is that migration is identical between all demes. This is likely to be problematic for species with limited vagility. In fact, for many species a model where migration occurs between neighbouring populations such as the stepping-stone is going to be more likely. At this stage it is unclear whether one could derive analytically the pdf of T 2 for a stepping-stone model. The work by [START_REF] Herbots | Stochastic models in population genetics: genealogy and genetic differentiation in structured populations[END_REF] suggests that it may be possible to compute it numerically by inversing the Laplace transform derived by this author. This work has not been done to our knowledge and would still need to be done.

Interestingly, this author has also shown that it is in principle possible to derive analytically the pdf of T 2 in the case of a two-island model with populations of different sizes. Again, this still needs to be done.

The SSPSC model has also been widely used [START_REF] Rogers | Population growth makes waves in the distribution of pairwise genetic differences[END_REF] and represents a first step towards using more complex models of population size change. For instance, the widely used method of [START_REF] Beaumont | Detecting population expansion and decline using microsatellites[END_REF] to detect, date and quantify population size changes [START_REF] Goossens | Genetic signature of anthropogenic population collapse in orangutans[END_REF][START_REF] Quéméré | Genetic data suggest a natural prehuman origin of open habitats in northern Madagascar and question the deforestation narrative in this region[END_REF][START_REF] Salmona | Signature of a pre-human population decline in the critically endangered Reunion Island endemic forest bird Coracina newtoni[END_REF] assumes either an exponential or a linear population size change. It should be straightforward to compute the pdf of T 2 under these two models because, as we explained above, the coalescent theory for populations with variable size has been very well studied [START_REF] Donnelly | Coalescents and genealogical structure under neutrality[END_REF]Tavaré 2004) and it is possible to write the pdf of T 2 for any demo-graphic history involving any type of population size changes. To go even further one could in principle use the history reconstructed by the PSMC [START_REF] Li | Inference of human population history from individual wholegenome sequences[END_REF] as the reference model of population size change and compare that particular demographic history to an n-island model using our results, and our general approach. At the same time, we should note that for complex models of population size change, including relatively simple ones such as the exponential model of [START_REF] Beaumont | Detecting population expansion and decline using microsatellites[END_REF], it is not straightforward to compute the number of differences between non recombining sequences. Significant work would probably be needed to apply the general approach outlined in our study to specific demographic models. But we believe that the possibilities opened by our study are rather wide and should provide our community with new interesting problems to solve in the next few years.

Comparison with previous work and generality our of results

The present work is part of a set of studies aimed at understanding how population structure can be mistaken for population size change and at determining whether studies identifying population size change are mistaken or valid [START_REF] Chikhi | The confounding effects of population structure, genetic diversity and the sampling scheme on the detection and quantification of population size changes[END_REF][START_REF] Heller | The confounding effect of population structure on Bayesian skyline plot inferences of demographic history[END_REF][START_REF] Paz-Vinas | The demographic history of populations experiencing asymmetric gene flow: combining simulated and empirical data[END_REF]. It is also part of a wider set of studies that have recognised in the last decade the importance of population structure as potential factor biasing inference of demographic [START_REF] Leblois | Genetics of recent habitat contraction and reduction in population size: does isolation by distance matter?[END_REF]Städler et al. 2009;[START_REF] Peter | Distinguishing between population bottleneck and population subdivision by a Bayesian model choice procedure[END_REF][START_REF] Chikhi | The confounding effects of population structure, genetic diversity and the sampling scheme on the detection and quantification of population size changes[END_REF][START_REF] Heller | The confounding effect of population structure on Bayesian skyline plot inferences of demographic history[END_REF][START_REF] Paz-Vinas | The demographic history of populations experiencing asymmetric gene flow: combining simulated and empirical data[END_REF] or selective processes [START_REF] Currat | Comment on "Ongoing adaptive evolution of ASPM, a brain size determinant in Homo sapiens" and "Microcephalin, a gene regulating brain size, continues to evolve adaptively in humans[END_REF]). Here we demonstrated that it is indeed possible to separate the SSPSC and StSI models. While we believe that it is an important result, we also want to stress that we should be cautious before extending these results to any set of models, particularly given that we only use the information from T 2 . Much work is still needed to devise new tests and estimation procedures for a wider set of demographic models and using more genomic information, including recombination patterns as in the PSMC method [START_REF] Li | Inference of human population history from individual wholegenome sequences[END_REF]. Beyond the general approach outlined here we would like to mention the study of [START_REF] Peter | Distinguishing between population bottleneck and population subdivision by a Bayesian model choice procedure[END_REF] who also managed to separate one structure and one PSC (Population Size Change) model. These authors used an ABC approach to separate a model of exponential PSC from a model of population structure similar to the StSI model. Their structured model differs from ours by the fact that it is not an equilibrium model. They assumed that the population was behaving like an n-island model in the recent past, until T generations in the past, but that before that time, the ancestral population from which all the demes derived was not structured. When T is very large their model is identical to the

StSI, but otherwise it may be quite different. The fact that they managed to separate the two models using an ABC approach is promising as it suggests that there is indeed information in the genetic data for models beyond those that we studied here. We can therefore expect that our approach may be applied to a wider set of models. We should also add that in their study these authors use a much larger sample size (25 diploid individuals corresponding to 50 genes). They used a maximum of 200 microsatellites which corresponds therefore to 10,000 genotypes, a number very close to the maximum number we used here. Altogether our study provides new results and opens up new avenues of research for the distribution of coalescent times under complex models.

Sampling and population expansions

Recent years have also seen an increasing recognition of the fact that the sampling scheme together with population structure may significantly influence demographic inference (Wakeley 1999) showed that in the n-island model genes sampled in different demes will exhibit a genealogical tree similar to that expected under a stationary Wright-Fisher model. Since our work was focused on T 2 we mostly presented our results under the assumption that the two genes of interest were sampled in the same deme. For diploids this is of course a most reasonable assumption. However, we should note that the results presented above also allow us to express the distribution of T 2 when the genes are sampled in different demes. We did not explore this issue further here, but it would be important to study the results under such conditions. Interestingly we find that if we assume that the two genes are sampled in two distinct demes, we can detect population expansions rather than bottlenecks. This could happen if we considered a diploid individual whose parents came from different demes. In that case, considering the two genes sampled in the deme where the individual was sampled would be similar to sampling his two parental genes in two different demes. Interestingly, this has also been described by

Conclusion: islands within individuals

To conclude, our results provide a general framework that can in theory be applied to whole families of models. We showed for the first time that genomic data from a single individual could be used to estimate parameters that have to our knowledge never been estimated. In particular we showed that we were able to estimate the number of islands (and the number of migrants) in the StSI model. This means that one can in principle use genomic data from non model organisms to determine how many islands make up the metapopulation from which one single individual was sampled. This is of course true as well for model organisms but it is particularly meaningful for species for which the number of individuals with genomic data is limited. Our ability to estimate n is one of the most powerful results of our study.

While such estimates should not be taken at face value, they surely should be obtained across species for comparative analyses. Also, during the last decade there has been a major effort to use programs such as STRUCTURE [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF] to estimate the number of "subpopulations" within a particular sample. Our work suggests that we might in principle provide additional results with only one individual. It is important to stress though that the answer provided here is very different from those obtained with STRUCTURE and similar methods and programs [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF][START_REF] Guillot | GENELAND: a computer package for landscape genetics[END_REF][START_REF] Chen | Bayesian clustering algorithms ascertaining spatial population structure: a new computer program and a comparison study[END_REF][START_REF] Corander | BAPS 2: enhanced possibilities for the analysis of genetic population structure[END_REF]). We do not aim at identifying the populations from which a set of individuals come. Rather we show that his/her genome informs us on the whole set of populations. In other words, even though we assume that there are n populations linked by gene flow, we show that each individual, is somehow a genomic patchwork from this (poorly) sampled metapopulation. We find these results reassuring, in an era where genomic data are used to confine individuals to one population and where division rather than connectivity is stressed.

  about here.] Structured population: Here we consider the classical symmetric n-island model (Wright 1931), see Figure 1 (b), where we have a set of n islands (or demes) of constant size N , interconnected by gene flow with a migration rate m, where M = N m is the number of immigrants (genes) in each island every generation. The whole metapopulation size is therefore nN (this is the total number of genes, not the effective size). Again, N is the number of haploid genomes, and N/2 the number of diploid individuals. We call this model the StSI, which stands for Structured Symmetrical Island model.

First moments :

 moments Equations 2 and 3 are different hence showing that it is in principle possible to identify genetic data produced under the two demographic models of interest. The two equations can be used to derive the expectation and variance of the two random variables of interest, T SSP SC 2 and T StSI 2

2

  and the T StSI 2 distributions, respectively (see Supplementary materials for the details of the estimation procedure). The ML estimates ( T , α) and ( M , n) were then used to define T SSP SC 2 or T StSI 2 reference distributions. The Kolmogorov-Smirnov (KS) test which allows to compare a sample with a reference distribution was then used to determine whether the observed n L sample could have been generated by the respective demographic models. In other words this allowed us to reject (or not) the hypothesis that the (t 1 , ..., t n L ) sample was a realization of the reference distributions (T StSI 2 or T SSP SC 2

Figure 4

 4 Figure4shows, for various values of n L , the results of the estimation of α (panels (a), (c), and (e), for simulations assuming α = 10 and T = (0.1, 1, 2), respectively ; see Supplementary Material for the other values) and the estimation of n (panels (b), (d), and (f) for simulations with n = 10 and M = (10, 1, 0.5), respectively; see Supplementary Material for the other values, corresponding to 26 figures and 168 panels). The first thing to notice is that both α and n are increasingly well estimated as n L increases. This is what we expect since n L

Figure 5

 5 Figure 5 shows for various values of n L the results of the estimation of T (panels (a), (c), and (e), for simulations assuming T = 0.2 and α = (2, 20, 100), respectively; see Supplementary Material for the other values) and the estimation of M (panels (b), (d), and (f), for

Figure 7

 7 Figure 7 is similar to Figure 6 but the data were simulated under the StSI model and the KS test was performed first using the SSPSC model as a reference ((a), (c), (e)) and then using the StSI model as a reference ((b), (d), (f)). The left-hand panels ((a), (c), and(e)) show results when we ask whether we can reject the hypothesis that genetic data were generated under a population size change model when they were actually generated under a model of population structure. In the right-hand panels ((b), (d), and (f)) we computed how often we rejected the hypothesis that genetic data were generated under the StSI model when they were indeed generated under that model of population structure. Altogether, the left-hand panels suggest that the results are generally best when M = (0.1, 0.2, 1), but that we get very good results for most values of M when we have 10, 000 loci and can reject the SSPSC when they were actually generated under the StSI model. The right-hand panels

Figure 8

 8 Figure 8 is divided in four panels showing the relationships between T and M (panels (b) and (d), for various values of α and n) and between α and n (panels (a) and (c),for various values of T and M ). In each of the panels we simulated data under a model for specific parameter values represented on the x-axis, and estimated parameters from the other model, and represented the estimated value on the y-axis. Since we were interested in the relationship between parameters (not in the quality of the estimation, see above), we used the largest n L value and plotted the average of 100 independent estimation procedures. In panel (a) we simulated a population size change (SSPSC) for various T values (represented each by a different symbol) and several values of α on the x-axis. We then plotted the estimated value of n for each case (i.e. when we assume that the data were generated under the StSI model).We find a striking linear relationship between these two parameters conditional on a fixed T value. For instance, a population bottleneck by a factor 50 that happened N 0 generations ago (T = 1) is equivalent to a structured population with n ≈ 22 islands (and M ≈ 0.71). Panel (c) is similar and shows how data simulated under a structured population generates specific parameters of population bottlenecks.Panels (b) and (d) show the relationship between T and M . We have plotted as a reference the curve corresponding to y = 1/x. As noted above and shown on this graph, this relationship is only approximate and depends on the value of α and n. Altogether, this figure exhibits the profound relationships between the model parameters. They show that the qualitative relationships between α and n, and between T and 1/M discussed above are indeed real but only correct up to a correcting factor. Still this allows us to identify profound relationships between population structure and population size change.

Figure 1 :

 1 Figure 1: Demographic models. (a): Single step population size change (SSPSC) model. The x-axis represents t, the time to the past in units of generations scaled by the number of genes. At time t = T , the population size changes instantaneously from N 1 to N 0 by a factor α. The y-axis represents the population sizes in units of N 0 (i.e. N (t)/N (0)). (b): Structured symmetrical island (StSI) model for n = 5 islands. Each circle represents a deme of size N . All demes are connected to each other by symmetrical gene flow, represented by the edges. In this example the total number of genes is 5N .

Figure 2 :

 2 Figure 2: Expected value and Variance of T 2 under the SSPSC and StSI models. This figure ilustrates how both models can have the same pair of values (E(T 2 ), V ar(T 2 )) for many sets of parameters. For the SSPSC model the time at which the population size change occured was fixed to T = 0.1 whereas α varied from 1 to 100 in one case, and T = 1.09, whereas α varied from 1 to 200 in the other case. For the StSI model the migration rate was fixed to M = 9 and M = 0.5, whereas n varies from 2 to 100.

Figure 3 :==

 3 Figure 3: Density of T 2 under the SSPSC and StSI models. Two sets of parameter values (panels (a) and (b), respectively) were chosen on the basis that expectations and variances were close. Panel (a): Density for the SSPSC model with T = 0.1 and α = 10.94, and for the StSI model with M = 9 and n = 10. For this set of parameters we have E T SSP SC 2 = 9.994, and E T StSI2

Figure 4 :

 4 Figure 4: Estimation of α and n. Panels (a), (c) and (e): Estimation of α under the SSPSC model for different sample sizes and T values. Simulations performed with α = 10 and T = (0.1, 1, 2). Panels (b), (d) and (f): Estimation of n under the StSI model for different sample sizes and M values. Simulations performed with n = 10 and M = (10, 1, 0.5).

Figure 6 :

 6 Figure 6: Proportion of rejected data sets simulated under the SSPSC model. Panels (a), (c) and (e): the reference model is the StSI model. Panels (b), (d), and (f): the reference model is the SSPSC, i.e. the model under which the data were simulated. Note that for the abscissa we used 2n L instead of n L because in order to perform the KS test it is necessary to first estimate the parameters using n L loci and then an independent set of n L values of T 2 .

Figure 7 :

 7 Figure 7: Proportion of rejected data sets simulated under the StSI model. Panels (a), (c), and (e): the reference model is the SSPSC. Panels (b), (d), and (f): the reference model is the StSI model, i.e. the model under which the data were simulated. Note that for the abscissa we used 2n L instead of n L because in order to perform the KS test it is necessary to first estimate the parameters using n L loci and then an independent set of n L values of T 2 .

  [START_REF] Peter | Distinguishing between population bottleneck and population subdivision by a Bayesian model choice procedure[END_REF] who found that when the 25 individuals were sampled in different demes, they would detect population size expansions rather than bottlenecks. Our results are therefore in agreement with theirs.Similarly Heller et al. (2013) also found that some signals of population expansion could be detected under scattered sampling schemes.

ACKNOWLEDGEMENTS

We are grateful to S. Boitard and S. Grusea for numerous and fruitful discussions and input throughout the development of this work. We thank I. Paz, I. Pais, J. Salmona, J. Chave for discussion and helpful comments that improved the clarity of the text. We thank J. Howard and the Fundação Calouste Gulbenkian for their continuous support. We are grateful to the genotoul bioinformatics platform Toulouse Midi-Pyrenees for providing computing and storage resources. This work was partly performed using HPC resources from CALMIP (Grant 2012 -projects 43 and 44) from Toulouse, France. This study was funded by the Fundação para a Ciência e Tecnologia (ref. PTDC/BIA-BIC/4476/2012), the Projets Exploratoires Pluridisciplinaires (PEPS 2012 Bio-Maths-Info) project, the LABEX entitled TULIP (ANR-10-LABX-41) as well as the Pôle de Recherche et d'Enseignement Supérieur (PRES) and the Région Midi-Pyrénées, France.

DISCUSSION

In this study we have analysed the distribution of coalescence times under two simple demographic models. We have shown that even though these demographic models are strikingly different (Figure 1) there is always a way to find parameter values for which both models will have the same first two moments (Figure 2). We have also shown that there are intrinsic relationships between the parameters of the two models (Figure 8). However, and this is a crucial point, we showed that the distributions were different and could therefore be distinguished. Using these distributions we developed a M L estimation procedure for the parameters of both models ( T , α) and ( M , n) and showed that the estimates are accurate, given enough genetic markers. Finally, we showed that by applying a simple KS test we were able to identify the model under which specific data sets were generated. In other words, we were able to determine whether a bottleneck signal detected in a particular data set could actually be caused by population structure using genetic data from a single individual. The fact that a single individual provides enough information to estimate demographic parameters is in itself striking (see in particular the landmark paper by [START_REF] Li | Inference of human population history from individual wholegenome sequences[END_REF]), but the fact that one individual (or rather sometimes as few as 500 or 1000 loci from that one individual) potentially provides us with the ability to identify the best of two (or more) models is even more remarkable. The PSMC (pairwise sequentially markovian coalescent) method developed by [START_REF] Li | Inference of human population history from individual wholegenome sequences[END_REF] reconstructs a theoretical demographic history characterized by population size changes, assuming a single non structured population. Our study goes further and shows that it could be possible to test whether the signal identified by the PSMC is due to actual population size changes or to population structure. However, models putting together these two scenarios have being proposed. In Wakeley (1999), a model considering a structured population who went through a bottleneck in the past is developed, allowing to estimate the time and the ratio of the bottleneck. Moreover, Wakeley (1999) discussed the idea that, in structured populations, changes in the migration rate and/or the number of islands (demes) can shape the observed data in the same way that impact of sampling schemes on the site frequency spectrum in nonequilibrium subdivided


Sjoerd Op 't Land
email: sjoerd.optland@eseo.fr

Senior Member, IEEE Mohamed ⇤ Ramdani
email: mohamed.ramdani@eseo.fr

Richard Perdriau
email: richard.perdriau@eseo.fr

Senior Member, IEEE Yannis Braux
email: yannis.braux@cst.com

M'hamed Drissi
email: mhamed.drissi@insa-rennes.fr

Using a Modified Taylor Cell to Validate Simulation and Measurement of Field-to-Shorted-Trace Coupling

Keywords: trace coupling, modified Taylor model, closed-form solution, full-wave simulation, GTEM cell

Predicting the immunity of electronic boards to radiated electromagnetic interference requires the computation of the coupling efficiency of an electromagnetic field to PCB traces. In the case of complex PCBs, full-wave electromagnetic solvers are convenient, yet at the expense of simulation time. Therefore, this paper introduces the extension of a modified Taylor-based analytical model to the case of traces terminated at one end by a non-characteristic impedance. This model makes it possible to determine the far-field-to-trace coupling using only a sum of closed-form equations. When applied to a shorted, meandered PCB trace, it was found to be as precise as 2.2 dB average absolute error with respect to GTEM measurements, which demonstrates its relevance for immunity prediction. Moreover, the full-wave simulation of this case study was validated using the extended model and found to be as precise as 1.4 dB average absolute error.

I. INTRODUCTION

Electromagnetic compatibility (EMC) problems can be understood as a three-element chain: source-transfer pathreceiver [START_REF] Paul | Introduction to Electromagnetic Compatibility[END_REF]. In the case of unshielded, wireless electronics, the dominant transfer path can consist in PCB traces. In that case, the routing of PCB traces is decisive for product compliance. Therefore, quickly estimating the field-to-trace coupling may be interesting.

Full-wave electromagnetic solvers can calculate this coupling for arbitrary trace geometry, terminal impedances, frequency, polarization and angle of incidence. However, these parameters need to be swept to find the average or worst case, which is time-consuming.

Therefore, other simplifications for this multi-dimensional problem have been discussed in the literature. Lagos developed an algorithm to find the worst case incidence for a straight trace with arbitrary terminal impedances and one frequency [START_REF] Lagos | Worst-case induced disturbances in digital and analog interchip interconnects by an external electromagnetic plane wave -part I: modeling and algorithm[END_REF]. Magdowski found closed-form equations for the average coupling under random illumination to a straight two-wire transmission line with arbitrary terminal impedances [START_REF] Magdowski | Closed-form formulas for the stochastic electromagnetic field coupling to a transmission line with arbitrary loads[END_REF]. Leone noticed that the coupling to matched terminal impedances is a reasonable approximation for moderate mismatches in terminal impedances [START_REF] Leone | On the coupling of an external electromagnetic field to a printed circuit board trace[END_REF].

Recently, a simple model was developed by the authors for a grazing, horizontally polarized incidence on a trace with matched terminal impedances: the modified Taylor cell [START_REF] Op 't Land | Simple, Taylor-based worst-case model for field-to-line coupling[END_REF], [START_REF] Op 't Land | Comparison of field-to-line coupling models: coupled transmission lines model versus single-cell corrected Taylor model[END_REF], which has a closed-form solution.

This paper introduces an improvement over this model for one mismatched terminal impedance. This model is then crosscompared with full-wave simulations and GTEM measurements on a realistic case study.

The paper is organized as follows: the aforementioned case study is first introduced in Section II: a meandered, shorted PCB trace illuminated by a grazing plane wave. The modified Taylor cell is then recalled in Section III. This model is then extended to include one non-characteristic load and is applied to the case study in Section IV. For reference, a fullwave electromagnetic simulation is set-up and applied to the case study in Section V. Finally, the experiment is performed in Section VI and compared with the full-wave solver and modified Taylor predictions. Conclusions are drawn and recommendations for future research are given in Section VII.

II. CASE STUDY:AMEANDERED,SHORTED MICROSTRIP

To allow for a meaningful comparison, a case study that can both be predicted theoretically and measured practically will be defined.

The theoretical prediction will be carried out both with modified Taylor cells and with a full-wave solver. As anything can be modeled in a full-wave solver, only the validity constraints of modified Taylor cells need to be considered. The modified Taylor cell predicts the coupling of a grazing, horizontally polarized plane wave to one straight, characteristic microstrip segment [START_REF] Op 't Land | Simple, Taylor-based worst-case model for field-to-line coupling[END_REF]. If the impedance discontinuities between segments are negligible, a finite cascade of microstrip segments of different planar orientation can be modeled by cascading Taylor cells [START_REF] Op 't Land | Comparison of field-to-line coupling models: coupled transmission lines model versus single-cell corrected Taylor model[END_REF].

Practically, the GTEM-cell was chosen as the plane wave source. The GTEM-cell is basically a tapered rectangular, paraxial 50 Ω waveguide with a metallic strip (the septum) as the center conductor, as schematized in Figure 1a. It sports a top aperture to illuminate a 10 ⇥ 10 cm PCB. The waveguide has a coaxially connectorized input, and its other end is terminated with a DC-18 GHz load, thus minimizing standing waves in the waveguide. As a result, the field strength under the PCB can be approximated by dividing the septum voltage V septum by the septum distance d [START_REF]integrated circuits -measurement of electromagnetic immunity 150 kHz to 1 GHz -part 2: measurement of radiated immunity -TEM-cell method[END_REF].

To measure the broadband transfer from the GTEM-cell input to any trace terminal, a Vector Network Analyzer (VNA) is used. Consequently, any monitored terminal is necessarily loaded by 50 Ω. The resulting measurement set-up is shown in Figure 1.

The PCB needs to provide coaxial connectors to give access to both trace terminals. Given the available surface, a threesegment microstrip is chosen for simplicity's sake. Because the far end has the stronger high-frequency coupling for a straight microstrip [START_REF] Op 't Land | Simple, Taylor-based worst-case model for field-to-line coupling[END_REF], the far-end induced voltage is studied. At the near end, a short calibration standard is connected, inevitably including an electrical delay ⌧ d . The entire path from connector to connector needs to be matched to the VNA characteristic impedance of 50 Ω. The resulting PCB and illumination are outlined in Figure 2.

III. MODIFIED TAYLOR MODEL

The analytical model used in this paper is based on Taylor's. This section recalls his model and the modification to take into account long line effects.

A. Taylor's Model

Taylor's model for field-to-line coupling represents the electric and magnetic induction as distributed sources along the transmission line. The specialization of the model for a twowire transmission line is shown in Figure 3. The component H n of the magnetic field (normal to the plane of the wires) induces a voltage, and the component E t of the electric field (in the plane and transversal to the wires) induces a current. Like a passive transmission line, each slice ∆z introduces a phase shift β∆z, where β is the wavenumber or spatial frequency.

If frequency asymptotically tends to zero, the exciting wavelength is infinitely long. Since the line has finite length, in this case, the phase shift along the line tends to zero and the field can be considered spatially uniform. Consequently, the line can be lumped as a single cell ∆z = `. In the case of characteristic loads (R ne = R fe = Z c ) the low-frequency eitherend terminal voltages can be found by inspecting Figure 3b [1]:

V LF = - 1 2 j! cE t Z c h`⌥ 1 2 j!µ 0 H n h`, ( 1
)
where c is the per-unit-length (pul) capacitance of the line. Unless otherwise noted, the near-end and far-end results are presented simultaneously throughout the paper; ⌥ means minus for the near end and plus for the far end.

B. Modification for Long Lines

If the field cannot be considered uniform along the line, the line needs to be lumped in infinitesimal slices dz, along each of which the field can be considered uniform. The terminal voltages can then be found by summing the contributions of this infinity of cells. Consequently, the resulting equations are no longer in closed form, which is inconvenient.

Fortunately, it turns out that the closed-form single-cell model of (1) can be maintained with a simple modification. A simple demonstration will be given here; the reader can refer to [START_REF] Op 't Land | Simple, Taylor-based worst-case model for field-to-line coupling[END_REF] for the formal derivation.

The starting point of this demonstration is a thought experiment on the lossless, characteristically terminated line of

E i ! c 0 H i 360 µm near end far end 1 2 3 `2 = 1.7 cm `1 = 2.5 cm `3 = 2.5 cm 50 Ω ⌧ d Figure 2.
Perspective on this paper's case study PCB: a three-segment, 50 Ω microstrip trace over a ground plane, illuminated by a horizontally polarized plane wave. The near end's coaxial connector is shorted, the far end's connector is terminated in a characteristic load (the VNA's input).

vertically Figure 3, illuminated from the near-end, that is, with the wave vector parallel to the line, pointing to the far end. The object of study is the far-end induced voltage, caused by a forward traveling wave on the line: the forward eigenwave. Its normalized amplitude w is just a phase lag:

w(z)=e -jβz ; β = ! v , (2)
where v is the phase speed of a wave on the transmission line.

Likewise, the illuminating field has a normalized amplitude i:

i(z)=e -jk p z ; k p = ! c 0 , (3)
where k p is the vector component of k parallel to the line.

Starting with an illumination frequency where the error introduced by only using a single Taylor cell is negligible, the following phenomena can be observed with increasing frequency along Figure 4(a-d).

Figure 4a shows the case for low frequencies, where the incident field remains the same along the line, so modeling the line as one Taylor cell is legitimate.

As shown in Figure 4b, the wavelength decreases. When the wavelength gets in the order of the line length, a propagating wave appears, both in free space and in the transmission line. Yet, this does not immediately invalidate the model. Indeed, the field is no longer uniform along the line, but the forward eigenwave of the line and the free space plane wave travel in the same direction. That means that, for every line slice, the free space wave and the eigenwave still have approximately the same phase. Therefore, it is still legitimate to model the line as one cell.

As frequency increases further, like in Figure 4c, the phase difference between the forward eigenwave and the incident plane wave becomes significant; in the example shown, the phase difference goes from 0 at z = 0 to ⇡ at z = `. On average, both waves are still cross-correlated, but less so than for low frequencies.

R ne R fe h

x z y (a) Line geometry: subscripts t, n and p (or x, y and z) denote excitation field components, respectively transversal, normal and parallel to the line segment. R ne and R fe are the near-end and far-end resistive terminations. In the extreme case of Figure 4d, the phase difference goes all the way from 0 at z = 0 to 2⇡ at z = `. On average, the two waves are no longer cross-correlated and coupling is not expected any longer.

+ - R fe R ne + - V ne + - j!µ0Hn h∆z j!cEt h∆z e -jβ∆z
To summarize, a single Taylor cell correctly predicts coupling for low frequencies. However, it overestimates the coupling at high frequencies. The single Taylor cell can be corrected for high frequencies with a measure for the lengthaverage cross-correlation between the line's eigenwave and the incident wave. This measure should be unity for low frequencies, in order not to modify the low frequency coupling. This unitless measure K should amount to zero when the phase difference along the line goes all the way from 0 to 2⇡.

The cross-correlation of the incident field and the line's eigenwave amplitudes is given by the complex conjugated product iw ⇤ . K is then found by averaging along the line:

K = 1 `Z 0 i(z) • w ⇤ (z) dz = 1 j(-k p ⌥ β)`⇣ e j(-k p ⌥β)`-1 ⌘ . (4)
To calculate the near-end induced voltage, the backward traveling eigenwave w = e +jβz is used. The resulting, modified Taylor cell is depicted in Figure 5.

IV. APPLICATION TO SHORTED AND MEANDERED TRACE

These modified Taylor cells will now be applied to predict the coupling to the shorted and meandered trace.

A. Near-End Non-Characteristic Load

In the case where the illuminated line is not terminated in a characteristic load at the near end, with a characteristic load still present at the far end, the reflected near-end voltage now adds up to the far-end induced voltage, as outlined in Figure 6a. Because of the far-end characteristic load, there are no further reflections to consider. This translates to Figure 6b and (5):

V fe | Γ ne 6 =0 = V fe | Γ ne =0 + Γ ne V ne | Γ ne =0 e -jβ`, (5)
where V fe is the far-end induced voltage and Γ ne is the voltage reflection coefficient (Z ne -Z c )/(Z ne + Z c ) of the near-end terminating impedance Z ne . For a single segment, the induced voltage can be found by inspecting Figure 5 with (4):

V fe | Γ ne 6 =0 = 1 2 j! (+µ 0 H n (0) -cE t (0)Z c ) h`K fe e -jβ`+ Γ ne 1 2 j! (-µ 0 H n (0) -cE t (0)Z c ) h`K ne e -jβ`. (6)
In the case study, the near-end impedance consists of a short circuit, preceded by a coaxial waveguide:

Γ ne = -1e -j!•2⌧ d , (7)
where ⌧ d is the one-way electrical delay of the coaxial connector and in the connected short-circuit standard.

B. Application to Meandered PCB Trace

To apply this modified Taylor cell to the case study, the two-wire line must be transformed into a microstrip trace [START_REF] Op 't Land | Simple, Taylor-based worst-case model for field-to-line coupling[END_REF] and a single segment into three segments [START_REF] Op 't Land | Comparison of field-to-line coupling models: coupled transmission lines model versus single-cell corrected Taylor model[END_REF].

The essential difference between a two-wire line and a microstrip is the presence of a ground plane and a substrate.

The ground plane doubles the field strength E i . This can be understood from the case of a Hertzian dipole, just above the zy-plane, infinitely far from the origin, such that the field E x is 1 V/m at the origin. If a ground plane is now placed at the zy-plane, there will be no field anymore under the zy-plane and the field above it will have doubled. It must be noted that this free space field 2E i equals the parallel plate field V septum /d of the GTEM-cell. As far as the dielectric substrate is concerned, the plane wave just above is imposed. Since the field in the substrate must follow with a constant phase lag, the wavenumber in the dielectric is equal to that in free space. Moreover, the material is not magnetic, hence the magnetic field H in the substrate is 2H i . Considering the substrate as part of an infinitely broad parallel-plate voltage divider, the electric field E in the substrate turns out to be 2E i /" r . These conclusions are summarized in Figure 7.

+ - Z c Z c + - V ne V fe + - e -jβj !µ 0 H n (0)h`K j!cE t (0)h`K
To model a three-segment microstrip, the contributions of the straight line segments `1, `2 and `3 must be summed:

V fe | Γ ne =0 =(K 1 fe V 1 LF,fe + K 2 fe V 2 LF,fe + K 3 fe V 3 LF,fe )e -jβ`( 8) V ne | Γ ne =0 =(K 1 ne V 1 LF,ne + K 2 ne V 2 LF,ne + K 3 ne V 3 LF,ne ), (9)
where K u fe is the far-end correction factor of the uth line segment. V u LF,fe (V u LF,ne ) denotes the the low-frequency far-end (near-end) voltage of the uth segment according to [START_REF] Paul | Introduction to Electromagnetic Compatibility[END_REF]. `is a shorthand for the line's total length `1 + `2 + `3.

The correction factors K u are found by integrating iw ⇤ along the line. Now that the line is no longer parallel with z everywhere anymore, let s denote the length along the line. For example,

K 2 = 1 `2 Z `1+`2 `1 i(s)w ⇤ (s) ds = 1 `2 Z `1+`2 `1
e -jk z `1 e ⌥jβs ds, (10) because the incident field is uniform along segment 2. In general,

K u = 1 j(-k s ⌥ β)`u ⇣ e j(-k z z end ⌥βs end ) -e j(-k z z begin ⌥βs begin ) ⌘ , (11)
where k s is the incident wave vector along the line segment and k z is the incident wave number.

Likewise, for the low-frequency contributions V u LF , the field orientation with respect to the line segment must be taken into account. For example,

µ 0 , ε 0 µ 0 , ε 0 ε r 2H i 2E i /ε r ! c 0 E i ! c 0 E i ! c 0 H i H i 2H i ! c 0 2E i
V 2 LF = -j!c E i " r Z c h`2 = -j ! c 0 E i h`2 p " r,eff " r , ( 12
)
because the incident magnetic field is parallel to line segment 2, hence it does not contribute to the terminal voltage.

Because the GTEM-cell has a 50 Ω input impedance and the far end is terminated in 50 Ω, the V fe /V septum ratio equals the transfer S 21 from septum to far end. A Python script took 3 ms on a Intel 2.53 GHz Core 2 Duo processor to calculate these far-end induced voltages according to [START_REF] Op 't Land | Simple, Taylor-based worst-case model for field-to-line coupling[END_REF]. The code can be viewed and run online [START_REF] Op 't Land | Extended Taylor online code and measurements for a meandered microstrip[END_REF], in comparison with the measurement data.

V. F ULL-WAV E SOLVER

The trace was entered in CST Microwave Studio as an infinitely-thin, perfect electric conductor (PEC). The ground plane was a 75 ⇥ 75 mm PEC beneath the trace. The substrate was a lossless dielectric layer with " r = 4.6, covering the ground plane. The trace was terminated in 50 Ω discrete face ports. The geometry was meshed using a hexahedral λ/10 mesh for 20 GHz, except near the microstrip: one mesh cell across the dielectric thickness and at least two mesh cells across the trace width. The Perfect Boundary Approximation (PBA) compensates for the staircase meshing of the mitered trace corners.

All boundary conditions were set to 'open and add space': that is, λ/8 extra space was added, terminated by 4 perfectly matched layers (PMLs). This geometry was then illuminated with a vertically polarized plane wave excitation, as illustrated in Figure 8a. The plane wave electric field strength was 1/d, where d is the GTEM-cell septum distance.

In order to take the terminal impedances into account, a geometry-circuit co-simulation was chosen, also because of its scaleability to many and non-linear loads. Hence, the above mentioned geometry with plane wave excitation was entered as one schematic block in the co-simulation. At the near end, a lossless transmission line was connected, representing the coaxial connection with delay ⌧ d up to the short standard, followed by a perfect short circuit. At the far end, a perfect load was connected, monitored by a voltage probe. As the structure is illuminated with the field caused by 1 V at the septum, the output voltage in volts predicts the transfer S 21 from septum to far end. The resulting schematic is depicted in Figure 8b.

The inherent causality of co-simulation requires a timedomain solver and a Finite Integration Technique (FIT)-based solver [START_REF] Weiland | A discretization model for the solution of Maxwell's equations for six-component fields[END_REF] was chosen. The co-simulation was then run with a -60 dB accuracy criterion on a 2 ⇥ 2.5 GHz Intel R E5 with 16 GB RAM and an NVIDIA R Tesla TM 20 graphics processing unit (GPU), which took 35 s. In order to allow for phase comparison, the phase reference of the plane wave had to be moved from the simulation boundary to the near end. Therefore, the simulated far-end voltage was advanced by the corresponding time of flight in the comparison script, of which the results are compared with the other methods in Figure 11.

The log-frequency average difference between simulation and the extended Taylor prediction, a measure for bias, was +0.8 dB over the 50 MHz to 20 GHz range. The average absolute error was 1.4 dB. Decreasing the accuracy criterion to -80 dB barely changed the simulation results, so the -60 dB accuracy criterion was kept. Decreasing the mesh size to λ/20 slightly improved the results: +0.4 dB average error and 0.9 dB average absolute error. Because this longer simulation yielded little improvement, the first simulation will be used in all comparisons.

VI. EXPERIMENT

This section describes the measurement set-up in more detail, discusses its quality and describes the measurement itself.

A. Set-up Details

In order to fabricate the case study with small and known mechanical uncertainties, the PCB was made in a four-layer industrial FR4 process, with mitered bends, traces matched to 50 Ω and connector-to-board transitions optimized, as can be seen in Figure 9. The PCB's ground contact was matched to the shape of the GTEM-cell aperture to smooth the transition from the ground plane to the GTEM-cell wall. The resulting PCB is shown in Figure 10. In order to express the layout as an objective function of fabrication and experiment constraints, it was described in code with PyPCB [START_REF]PyPCB -specifying PCB layout in Python[END_REF].

To numerically predict the measurements, some set-up specific numbers had to be entered into both the analytical model and the full-wave simulation. The GTEM septum distance d was measured to be 42.2 mm on average under the PCB aperture. From the near end of the trace to the short circuit, the signal traverses the PCB, the SMA connector and the short standard. The PCB has a thickness t of 1.55 mm and a relative permittivity " r of 4.6. Using Agilent's ADS, the filling factor of the microstrip was calculated to be 0.67 on average over the frequency range of interest. With this approximation, the microstrip's effective relative permittivity " r,eff was calculated to be 3.4. Modeling this path as a homogeneous waveguide, the delay p " r • t/c 0 should approximate 11.1 ps. The SMA connector has a Teflon dielectric (" r = 2.1) of 7.67 mm from board to reference plane. Similarly, this should account for 37.1 ps. The Agilent 85052D-60006 short circuit is specified to have a 31.78 ps delay, totaling to a delay ⌧ d of 79.95 ps.

B. GTEM Field Uniformity

Before illuminating the PCB with the GTEM cell, the field uniformity needs to be assessed. Therefore, E-field and H- 5). Port 1 is the GTEM cell input, port 2 is the trace's far end, while the near-end is terminated in a short circuit. This figure can be reproduced online [START_REF] Op 't Land | Extended Taylor online code and measurements for a meandered microstrip[END_REF].

field sensors were integrated in one corner of the PCB. The H-field sensor consists of two vias and a trace that form a loop tangential to the septum, as short as technologically possible for the field to be uniform in the loop (1.23 mm), as high as technologically possible to increase its surface (1.50 mm). To measure the E-field, a disk of 3 mm diameter was used, in order to stay under λ/10 at 10 GHz and to have sufficient coupling surface to measure strong enough a signal.

The transfer from GTEM cell septum to the E-field sensor was measured at four different PCB orientations by rotating the PCB. Up to 5 GHz, the log-frequency averaged difference between the maximum and minimum transfer was 1.4 dB, and 5.8 dB between 5-20 GHz. As far as the H-field sensor is concerned, the average difference was 1.5 dB up to 5 GHz, and 3.9 dB between 5-20 GHz, taking into account both positions. These numbers suggest that the field non-uniformity is in the order of a few dB.

C. Measurements

The measurement of Figure 1 was performed, yielding the S 21 plotted in Figure 11. The log-frequency average difference between measurement and extended Taylor prediction was +0.2 dB over the 50 MHz to 20 GHz range. The average absolute error was 2.2 dB.

The ripple in the measurement is supposedly caused by standing waves in the GTEM-cell. However, no satisfactory explanation was found for the low-frequency offset of +2 dB with respect to both the modified Taylor model and the fullwave simulation; explanations in the literature deal with zerocentered deviations from the V septum /d estimator. In general, it should be noted that quite some imperfections of the measurement set-up were not taken into account in model and simulation: the substrate permittivity being frequencydependent, skin loss and dielectric losses in the traces, the discontinuity between the GTEM-cell wall and the substrate, and the non-uniformity of the illuminating field.

VII. CONCLUSIONS AND PERSPECTIVES

The coupling of an incident grazing plane wave to the farend of a meandered, near-end shorted PCB trace with ground plane was studied. To that end, the modified Taylor model was extended to include a near-end non-characteristic load. This extension should hold for any electrically long, lossless transmission line, illuminated by a linearly polarized plane wave and terminated characteristically at its far end. The solution was still closed-form, so the calculation time was only 3 ms.

Then, a full-wave electromagnetic solver was compared against this model. In CST, an infinitely-thin, perfectlyconducting microstrip trace with a lossless dielectric was illuminated with a far-field source. Using time-domain cosimulation, the near-end short circuit was taken into account. The solution differed by +0.8 dB from the extended Taylor prediction, averaged over log-frequency in the 50 MHz to 20 GHz band. The average absolute error was 1.4 dB. The simulation time was 35 s.

Finally, a measurement in a GTEM cell was performed and compared against the extended Taylor model. The average difference with the extended Taylor model was +0.2 dB with an average absolute error of 2.2 dB. This error is in the order of the field strength variation, as measured with integrated H-and E-field sensors. As the correlation with measurement was already quite good, dielectric losses, amongst others, were not incorporated in the full-wave simulation neither in the modified Taylor model.

To sum up, the extended Taylor model, despite its many simplifying assumptions, is able to quickly estimate the field-totrace coupling for grazing incidence and one non-characteristic load.

The next step would be to adapt the modified Taylor model to approach the near-field illuminated suspended wire of IEEE benchmark case 5.3.1.2 [START_REF][END_REF]. This requires to take into account a far-end non-characteristic load (1 kΩ), which would introduce the complexity of bidirectional reflections and require another measurement method. For industrial applicability, nongrazing, arbitrary polarized incidence should be taken into account as well, probably using worst-case analysis to simplify the solution. Finally, the solution should be made accessible to PCB designers, to estimate PCB immunity quickly and early in the design process.

2 (

 2 GTEM frontal crosssection with electrostatic (solid, arrows) and magnetostatic (dashed) field lines. VNA 1 b) Overview lateral cross-section with electrostatic field lines. The resistive/absorbent hybrid load provides broadband characteristic termination. (c) Overview photo, relevant parts highlighted.

Figure 1 .

 1 Figure 1. Measurement set-up: GTEM-cell, VNA and test PCB.

  V fe (b) Taylor's cell: approximation of passive transmission line slice ∆z, with a voltage source representing the electromotive force (emf) and current source representing the electrostatic force (electric induction). c denotes the per-unitlength (PUL) capacity of the line, β is the wavenumber of the line.

Figure 3 .

 3 Figure 3. Modeling the coupling of an electromagnetic wave to a two-wire transmission line.

Figure 4 .

 4 Figure 4. Phase along the transmission line of the line's eigenwave ∠w (solid line) and illuminating plane wave ∠i (dashed line), for increasing frequency.

Figure 5 .

 5 Figure 5. Modified Taylor cell taking into account long-line effects.

Figure 6 .

 6 Figure 6. Propagating waves and induced voltage calculation for the case of a non-characteristic near-end load and a characteristic far-end load.

Figure 7 .

 7 Figure 7. The far-away plane wave source (top left) is reflected by the microstrip's ground plane (image source at the bottom left). This results in the fields shown at the right.

  Geometry (figure not to scale, in order to see the substrate) as entered into CST Microwave Studio. The dielectric substrate is shown translucent. Note the polarization of the incident plane wave from the left. (b) Co-simulation schematic as entered into CST Design Studio. The first block is a lossless transmission line defined by its delay ⌧ d , the second block is the geometry with excitation of Figure 8a.

Figure 8 .

 8 Figure 8. Overview of the CST Studio full-wave simulation.

  Frontal cross-section. The M2 screws are made of nylon.

Figure 9 .Figure 10 .

 910 Figure 9. SMA-to-microstrip launch: a Molex SD-73251-185 SMA jack, screwed onto a footprint designed using CST to minimize reflections (drawing not to scale, dimensions in mm).

Figure 11 .

 11 Figure 11. Field-to-trace coupling of the meandered trace, illuminated in a GTEM cell, as measured, full-wave simulated and calculated with the extended Taylor model (Equation5). Port 1 is the GTEM cell input, port 2 is the trace's far end, while the near-end is terminated in a short circuit. This figure can be reproduced online[START_REF] Op 't Land | Extended Taylor online code and measurements for a meandered microstrip[END_REF].

ACKNOWLEDGEMENTS

This research project was co-financed by the French national project SEISME (simulation of emissions and immunity of electronic systems). We thank Mohamed Latrach for his effort to provide us with state-of-the-art tools and equipment.


Anne Jost

Médard Thiry

Agnès Rivière

Sophie Violette

Julio Gonçalvès

Groundwater silicification as a proxy of paleo-permafrost depth and a constraint for a fluid flow and geothermal modelling

Introduction

In the Paris basin, the Fontainebleau Sands (Early Oligocene) contain superposed flat-lying lenses of very tightly cemented sandstones, the origin of which was related to paleo-groundwater (Thiry et al. 1988). Recent dating of calcite crystallarias, sometimes included in the quartzites, has provided an age constraint that refers to cold periods of the Quaternary. It suggests a new mechanism of sand cementation, favoured by the decrease in silica solubility when feeding warm groundwater would discharge close to a freezing front (Thiry et al. 2013). Quartzite lenses would therefore act as a proxy for the presence and thickness of permafrost. Model calculations aims to investigate the temperature and pressure conditions associated with a cold front propagation near the groundwater outlet. We discuss them in terms of constraints on fluid flow and geothermal profile to provide sustainable conditions to achieve silica deposition.

Mechanisms of silicification

The discontinuous distribution of the silicified bodies, as well as the correlation between their localisation and the recent or present morphology, suggest a relatively recent surficial silicification, near the outcrop zones.

The general arrangement of the sandstone in subhorizontal layers and their elongated morphologies towards valleys may also indicate a control on their genesis by paleo-groundwater (Thiry et al. 1988, fig 2). Each sandstone level would correspond to an old groundwater level and form at the groundwater discharge zones. Thiry et al. (1988) thus associated the successive sandstones development in the sand formation to the falling water table levels during the downcutting of the valleys, in a first silicification model. ( )

dt dS ρ ρ t p g S g ρ z p μ T kk ρ z x p μ T kk ρ x i i w s w r w r w - + ∂ ∂ =         + ∂ ∂ ∂ ∂ +         ∂ ∂ ∂ ∂ φ ) ( ) ( ( ) ( ) t S L ρ t T ρC Tq C ρ T i i w w e ∂ ∂ - ∂ ∂ = ⋅ ∇ - ∇ ⋅ ∇ φ κ Loing Riv.

Corbeil

Discussion & Conclusions

The duration of silicification is geologically short but in view of the weak solubility of silica in surficial waters, substantial groundwater flow is nonetheless needed to supply the silica precipitated from the solution. It shows the importance of a high groundwater flow rate to sustain silica for sandstones cementation in a time span compatible with the hydrothermal constraints. Near the groundwater oultets, where present springs flow at a rate of 3-150 m 3 /hr, higher velocities would allow a convective thermal regime to maintain the warmer temperatures of the sand aquifer and a sharp geothermal gradient on contact with discontinuous permafrost, expected to provide favourable conditions to achieve silica deposition. An abrupt transition towards a cold period is somehow also needed (our modelling assumes it is instantaneous) to prevent the geothermal gradient from slowly adjusting to decreasing surface temperatures. However, the presence of a large vadose zone on the plateaus (not yet taken into account in the numerical experiments) may play an additional insulator role.

Future research will be conducted accordingly. Fluid and thermal conditions during thawing of permafrost should also be considered.

Running a geochemical model would eventually allow to refine our hypothesis and verify if several hundreds of years of silica precipitation in these conditions are sufficient to generate the Fontainebleau Sandstones.

  Fig. 1 -Hydrogeologic framework (a) Geological map of the Beauce Plateau with location of the cross-section line. (b-c) Typical 2-D hydrogeologic section of the Beauce Plateau with its simplified version used for modelling. Silicification in the Fontainebleau Sands produced up to 4 principal levels of superposed flat-lying lenses of very tightly cemented sandstones (0.5 to 8 m thick, m to hm long, see fig 2). Drill-hole data indicate that these sedimentary quartzite layers are restricted to outcrops on the valley slopes and do not extend more than a few hundred meters beneath the overlying limestone cover of the plateaus (fig 1b). The sandstone lenses are discontinuous and of variable size, from a few dm 3 volume to several thousands of m 3 .

Fig. 2 -

 2 Fig. 2 -Sandstones, Darvault (77) Contorted morphologies of the sandtone lenses with often asymmetrical profiles, elongated towards valleys and interpreted as evidence of the paleo-groundwater flow direction.

Fig

  Fig. 3 -Regional distribution of the Sandstones

Fig. 4 -

 4 Fig. 4 -Preferential cementation of sands along dune ridges resulting in parallel stripes of sandstone layers

Fig

  Fig. 5 -Calcite crystallarias, Fontainebleau Sands, Larchant (77)

Fig. 7 -

 7 Fig. 7 -Model of silicification associated with a cold front propagation During cold periods of the Quaternary, silicification would occur at the abrupt temperature interface between permafrost and warmer groundwater.

  -NNE cross-section (70 km long, 125 m thick, fig 1c) of the Beauce Plateau serves as the basis for 2-D coupled fluid flow and heat transport simulations. The data used are from various studies on the Fontainebleau Sands and the Beauce aquifer: We used a newly developed numerical code that simulates groundwater flow and heat transfer in the (x,z) plane with integrated freezing and thawing processes (latent heat effects and modifications of hydraulic and thermal properties due to ice formation; Rivière et al. 2013). It solves the following equations: Groundwater flow Heat transport Steady state simulations were run under modern conditions (annual surface temperature of 11°C, recharge of ~10 mm/yr) to generate initial conditions, characterized by a vadose zone up to several tens of meters and a groundwater flow rate of 1-10 m/yr. A continental divide crosses the section, between the Loire and Seine catchments (fig 8a). An average geothermal gradient is calculated as 30.7°C/km in the Beauce plateau for a uniforme basal heat flow of 70 mW.m -2 (see fig 9). With a modified Peclet number Pe* at the aquifer scale (van der Kamp 1984) widely lower than 0.1, groundwater flow does not affect the thermal regime (fig 8a). How long a relatively warm groundwater flow could have prevailed in the Beauce aquifer? Starting from modern conditions, the top boundary condition was altered to simulate a cold surface temperature (-0.5°C). It is assumed that frozen ground prevents recharge from the surface; however positive groundwater outlet temperatures are maintained. As a result, hydraulic heads decrease slowly (fig 8b). A thin layer of frozen ground (1-2 m thick after 500 yr, fig 8b) developed on top of the plateau, spatially depending on the heterogeneous properties of sands and limestones. Geothermal profile (fig 9) varies accordingly. Temperature as high as 11°C can be maintained during at least ~500 kyr (fig 8b & 9) at the base of the sands aquifer (~60 m deep) upstream from the edges of the plateau.

Fig

  Fig. 9 -Geothermal profile evolution in the course of the simulation Location: ~6 km upwards from the main outlet; time in kyr.

  Initial temperature distribution and water table (dashed line).

Fig. 8 -

 8 Fig. 8 -Model results