Carlos A Pires
email: capires@fc.ul.pt

Olivier Talagrand

Marc Bocquet

Diagnosis and impacts of non-Gaussianity of innovations in data assimilation

Keywords: Non-Gaussianity, Data Assimilation, Innovation Statistics, Maximum Entropy Method PACS: 92.60.-e, 93.65.+e, 93.85.+e, 92.90.+x

The Best Linear Unbiased Estimator (BLUE) has widely been used in atmospheric and oceanic data assimilation (DA). However, when the errors from data (observations and background forecasts) have non-Gaussian probability density functions (pdfs), the BLUE differs from the absolute Minimum Variance Unbiased Estimator (MVUE), minimizing the mean square a posteriori error. The non-Gaussianity of errors can be due to the inherent statistical skewness and positiveness of some physical observables (e.g. moisture, chemical species) or because of the nonlinearity of the data assimilation models and observation operators acting on Gaussian errors. Non-Gaussianity of assimilated data errors can be justified from a priori hypotheses or inferred from statistical diagnostics of innovations (observation minus background). Following this rationale, we compute measures of innovation non-Gaussianity, namely its skewness and kurtosis, relating it to: a) the non-Gaussianity of the individual error themselves, b) the correlation between nonlinear functions of errors, and c) the heteroscedasticity of errors within diagnostic samples. Those relationships impose bounds for skewness and kurtosis of errors which are critically dependent on the error variances, thus leading to a necessary tuning of error variances in order to accomplish consistency with innovations. We evaluate the sub-optimality of the BLUE as compared to the MVUE, in terms of excess of error variance, under the presence of non-Gaussian errors. The error pdfs are obtained by the maximum entropy method constrained by error moments up to fourth order, from which the Bayesian probability density function and the MVUE are computed. The impact is higher for skewed extreme innovations and grows in average with the skewness of data errors, especially if those skewnesses have the same sign. Application has been performed to the quality-accepted ECMWF innovations of brightness temperatures of a set of High Resolution Infrared Sounder (HIRS) channels. In this context, the MVUE has led in some extreme cases to a potential reduction of 20-60% error variance as compared to the BLUE.

Introduction

Most of the operational atmospheric-oceanic data assimilation systems rely on linear estimation theory that aim at to compute the Best Linear Unbiased Estimator (BLUE) of the system state from all available data (observations and background) and the corresponding error covariance matrices. For some cases, errors of assimilated data (direct or indirect observations and a priori background estimates) have non-Gaussian probability density functions (pdfs). When data errors are non-Gaussian, the a posteriori error variance minimiser is the Minimum Variance Unbiased Estimator (MVUE), obtained from the mean of the Bayesian a posteriori pdf combining all data [START_REF] Talagrand | Bayesian estimation. Optimal interpolation. Statistical linear estimation[END_REF] . Nevertheless, the MVUE restricted to the subspace of linear estimators coincides with the BLUE, even when data errors are Non-Gaussian.

Sampling statistics and histograms of innovations (observation minus projected background ) often reveal statistically significant non-Gaussianity, thus proving that observation or background errors are separately non-Gaussian or even that errors are jointly non-Gaussian in the observational space. This in turn gives proof of the BLUE sub-optimality as far as the a posteriori error variance minimization is concerned.

The assumption of Gaussian errors generally comes with the linearity of the observation operators and Gaussianity of errors in state variables, which is not true in general. In fact, assimilated observations can have instrumental and representativeness errors of non-Gaussian nature, especially if they come from complex discontinuous and/or nonlinear observation operators such as in the cases of specific humidity, tracing and reactive chemical species (e.g. ozone) [START_REF] Errera | Four-dimensional variational chemical assimilation of crista stratospheric measurements[END_REF][START_REF] Elbern | Ozone episode analysis by four-dimensional variational chemistry data assimilation[END_REF] cloud droplet concentrations, land-species concentrations, remote sensing data (e.g. satellite radiances, scatterometer and Doppler radar data).

Moreover, Gaussianity of errors implies a symmetric probability density function (pdf) around the mean, which is not valid for example in cases dealing with moisture, gaseous pollutants and aerosols. There are also extreme events like droughts and floods, which are not accommodated by a Gaussian pdf. Another source of non-Gaussianity in observation errors is the presence of extreme gross errors [START_REF] Ingleby | Bayesian quality control using multivariate normal distributions[END_REF]. Data assimilation schemes coping with non-Gaussian error pdfs such as the lognormal are now available [START_REF] Fletcher | A hybrid normal and lognormal distribution for data assimilation[END_REF][START_REF] Fletcher | A data assimilation method for log-normally distributed observational errors[END_REF][START_REF] Bocquet | High resolution reconstruction of a tracer dispersion event: application to etex[END_REF]. When errors are non-Gaussian, one can still formulate a cost-function with the minimum verifying at the maximum of the a posteriori Bayesian pdf, also designated by the Maximum A Posteriori (MAP). The referred cost-function can be decomposed into a quadratic term corresponding to Gaussian errors, with the minimum at the BLUE, and an additional term measuring the departure from Gaussianity. The cost function minimum, the MAP, is shifted from the BLUE in accordance to the weight of non-Gaussianity. This appears as a correcting term due to a priori ill-specified error pdfs or equivalently to a generalized parametric bias. The removal of parametric bias, as considered above, is automated by the Variational Bias Correction (VarBC) [START_REF] Dee | Bias and data assimilation[END_REF] by the addition of

Description of the Data

The data on which the diagnostics are performed below is extracted from the database of ECMWF. It consists on quality-controlled observations of brightness temperature (BT) performed over the ocean by the High-resolution Infrared Sounder (HIRS) carried by satellite NOAA-17 (http://www.ncdc.noaa.gov). The observations have been performed in 8 wavelength channels: 4, 5, 6, 7, 11, 12, 14 and 15, over 8 different days in the period October 2005-January 2006: 20051015, 20051030, 20051115, 20051201, 20051215, 20051230, 20060115 and20060130 in yyyyMMdd notation. The main statistical features of the data errors are given in Table 1. In the ECMWF assimilation system, both the observations and the background are assumed to be unbiased. The specified standard deviation σ so of observational error is constant for each channel. The standard deviation σ sb of the background error as specified by the ECMWF assimilation scheme depends on the channel, pixel position and temporal instant [START_REF] Derber | A reformulation of the background error covariance in the ecmwf global data assimilation system[END_REF][START_REF] Anderson | Diagnosis of background errors for radiances and other observable quantities in the variational data assimilation scheme and the explanation of a case of poor convergence[END_REF]. The five channels: 5, 6, 7, 14 and 15, are sensitive to sea-ice, whose temperature is considered to be ill-predicted by the forecast model. The corresponding background error standard deviation is therefore taken to be larger in the conditions when sea-ice is forecast to be present, preferentially at latitudes poleward of 60 o . Fig. 1, relative to channel 6, shows this fact through the scatter plot of pairs (σ sb , North latitude), thus putting in evidence the clustering in the distribution of σ sb . This has led to partitioning the data pairs (observation, background) for each of referred above five channels into two subsamples, meant to correspond respectively to 'ice-free' and 'ice-covered' surface conditions. In the case of Fig. 1 for channel 6, the partitioning is defined by the thick lines; a similar partitioning is used for the other ice-sensitive channels. Samples corresponding to 'ice-free' surface conditions are identified by the subscript s: 4s, 5s, 6s, 7s, 11s, 12s, 14s and 15s as in Table reftable1, where the leading number stands for channel. Data points under 'ice-free' surface conditions were previously 'thinned' to a mutual distance of at least 300 km. The samples for 'ice-covered' conditions, identified by the subscript i: 5i, 6i, 7i, 14i and 15i were not subject to any thinning in order not to drastically reduce the sample size. For each of the 13 presented data samples, N denotes the size of the working samples while m(σ sb ) and std(σ sb ) denote, respectively, the mean and standard deviation of σ sb estimated from data. It is seen that both m(σ sb ) and std(σ sb ) are larger for the 'ice-covered' samples, comparatively to the 'ice-free' samples for each of the channels 5, 6, 7, 14 and 15.

Innovation Statistics

We recall that the innovation corresponding to an observation y = x + ε o , where x is the exact value and ε o , is the observational error, is defined as the difference:

d = y -x b = ε o -ε b , (1)
where x b = x + ε b is the background estimate for x, and ε b the corresponding error. The innovation is the only combination of the data (observation and background) that is independent of the exact value x, and is as such the only objective source of information on the errors ε o and ε b . However, it clearly brings no information on either ε o or ε b individually.

Bias and Variance

Rows 5 and 6 of Table 1 show the bias b d and the standard deviation σ d of the innovations, computed along each set. For all the 13 samples, the bias b d , is of the order of a few hundredths of a Kelvin degree, while the standard deviation σ d is typically larger by one or two orders of magnitude. The biases are statistically significant, as shown by the fact that the ratio |b d |/(σ d / √ N ) is significantly larger than 1; it varies between 10 and 40.

In the ECMWF assimilation system, background and observation errors are assumed, as is usual, to be mutually uncorrelated. The variance of the innovation should therefore be the sum of the variances of the two errors. This is not so in the present case. It can be seen from the Table 1 that the latter sum is larger than the innovation variance for all 13 subsets, by a large factor which varies between 1.84 (sample 11s) and as much as 27.72 (sample 7i). This is due in part to the assumed variance for observational error, which (with the only exception of sample 11s) is larger than the variance of the innovation. This corresponds to a purposeful inflation of the variance of the observational error in the ECMWF assimilation system. That inflation is meant to compensate for the fact that the spatial correlation of errors in satellite observations is ignored. The excess of total error variance is also due to the assumed variance for the background error in the 5 'ice-covered' samples, which is significantly larger than the variance of the innovation.

Higher-order Statistical Moments and non-Gaussianity

The classical skewness and kurtosis of innovations defined as:

s d ≡ (d /σ d ) 3 (2a)
k d ≡ (d /σ d ) 4 -3, (2b)
are respectively the third and fourth-order standardized cumulants of innovations, where the bar and the prime stand respectively for (channel+surface condition) sampling mean and departure with respect to the mean. The skewness can take all values in the interval]-∞, +∞[. When s d = 0, the corresponding pdf is asymmetric with respect to its mean, with in general a more extended tail of extreme values on the side with the same sign as s d . For a Gaussian pdf, s d = 0 and k d = 0, so that any statistically significant deviation from those values is a proof of non-Gaussianity. The probability distribution is said to be platykurtic (leptokurtic) if the kurtosis k d is negative (positive). It thus has tails that are less populated (more populated) than those of the Gaussian distribution with the same mean and variance. The kurtosis is very sensitive to the presence of outliers, which can significantly increase its value. There is a lower bound for k d as a function of s d [START_REF] Wu | Calculation of maximum entropy densities with application to income distribution[END_REF].

The kurtosis ranges in the interval [-2, +∞[. This is seen from the variable

r d ≡ (d /σ d ) 2 -s d (d /σ d ), which
has sampling variance:

η d ≡ var (r d ) = k d -s 2 d + 2 ≥ 0. ( 3
)
The The relative entropy of a pdf ρ 1 (x) with respect to another pdf ρ 2 (x), also called the Kullback-Leibler divergence between ρ 1 (x) and ρ 2 (x), is defined as:

KL(ρ 1 , ρ 2 ) ≡ ρ 1 (x) ln ρ 1 (x) ρ 2 (x) dx. ( 4
)
It is equal to zero for ρ 1 = ρ 2 and is strictly positive for ρ 1 = ρ 2 . It also possesses the fundamental property of being invariant in any one-to-one mapping x ↔ ξ. Despite not defining a distance in the sense of metric spaces (it is not symmetric), the relative entropy is a useful measure of the discrepancy between ρ 1 and ρ 2 .

The negentropy of a pdf ρ in the sense of information theory, is the relative entropy of ρ with respect to the Gaussian pdf ρ G with the same expectation and variance [START_REF] Cover | Elements of Information Theory[END_REF]. For small absolute values of skewness s d and kurtosis k d , the negentropy is approximated thanks to the Gram-Charlier expansion of ρ/ρ G [START_REF] Comon | Independent component analysis, a new concept?[END_REF]:

KL(ρ, ρ G ) ≈ 1 12 s 2 d + 1 48 k 2 d ≡ J d . (5)
The approximate negentropy J d , expressed in terms of s d and k d , will be used as an integrated diagnostic of non-Gaussianity.

Rows 7 to 9 of Table 1 show the three quantities s d , k d and J d for the 13 samples of innovations. For each sample, checks of deviation from Gaussianity have been performed on the basis of an ensemble of 1000 realizations, each with the same size N , of the Gaussian probability distribution with the same expectation and variance of the sample. Values marked in bold in rows 7 to 9 of the Table 1 indicate non-Gaussianity at the corresponding 99% level of confidence. It is seen that only sample 14s can be considered Gaussian by those checks and 10 subsets out of 13 are non-Gaussian by all three checks. It can also be noted that the kurtosis is positive (leptokurtic) for the 8 'ice-free' samples. The 5 'ice-covered' samples are platykurtic (negative kurtosis) and negatively skewed. In order to show a graphic evidence of non-Gaussianity, we present in Fig. 2, the histogram of innovations for the sample 4s which is the more skewed one or having the largest |s d |. In that Figure, we superpose the Gaussian pdf fit and the non-Gaussian pdf fit obtained by the ME method, constrained by moments up to the fourth order. There is a clear difference between the two fits. It is interesting to elucidate the origin of those features as result of the observations and background statistics, the quality control and other possible sources of non-Gaussianity. The elimination, through the quality control, of too high absolute values of innovations can lead to an extra reduction of innovation extremes and therefore leading to possible platykurtosis of innovations. It is seen that there are significant inconsistencies between the textita priori hypotheses that lie at the basis of the assimilation system on the one hand, and the observed statistics of the innovation of the other. Both the first and the second-order statistical moments of the innovation are inconsistent with the a priori specified biases and variances of the observation and background errors. In addition, the innovation turns out to be significantly non-Gaussian, meaning that the assimilation system, even if optimal in the sense of the BLUE, does not achieve Bayesian estimation. These conclusions are not limited to the satellite observations considered above. Similar results

(not shown) have been obtained in the ECMWF system with radiosonde observations. The only significant difference with the case of satellite observations is that the amplitude of the observational error does not seem to be overevaluated (or at least not as systematically).

Statistical consistency between innovations and errors

Mean and Variance

The statistical moments of the innovations d can be expressed as combinations of single and joint moments of the errors ε o , ε b , leading to consistency relationships between the statistics of errors, which can be specified either by the assimilation scheme or as additional error hypotheses. The statistical moments can also be expressed in terms of cumulants [START_REF] Kenney | Cumulants and the cumulant-generating function, additive property of cumulants, and sheppard's correction[END_REF][START_REF] Maurice | Kendall's Advanced Theory of Statistics[END_REF], thus leading to related consistency cumulant-relationships. In order to accomplish consistency, we need in general to perform some tuning to the specified error statistics leading in particular to bias error correction and tuning of error variances [START_REF] Desroziers | Diagnosis and adaptative tuning of information error parameters in variational assimilation[END_REF][START_REF] Chapnik | Properties and first application of an error-statistics tuning method in variational assimilation[END_REF]. Therefore, let us consider the statistically consistent biases b co and b cb , respectively of ε o and ε b , taken uniform over each sample. Then we introduce tuning factors (f o , f b ) taken constant within each innovation sample. These factors relate the consistent standard deviations of the unbiased errors ε o and ε b , noted (σ co , σ cb ), to the specified standard deviations of errors (σ so , σ sb ) as: σ co ≡ f o σ so and σ cb ≡ f b σ sb . The values of σ sb and σ cb vary within samples meaning that background errors as well as innovations own heteroscedasticity, whose presence has to be taken into account in the consistency relationships. In order to isolate its effect, we have considered the standardized observation and background errors:

ε * o ≡ ε o -b co σ co ≡ ε o σ co , (6a)
ε * b ≡ ε b -b cb σ cb ≡ ε b σ cb . (6b)
In this paper, we will make the fundamental hypothesis that the joint and marginal assumed pdfs of the standardized errors ε * o and ε * b , as well as the biases b co , b cb , the standard deviation of observation error σ co , and the tuning factors f o , f b , are independent both from the observation pixel i.e. its geographical position and from the observation time, being considered uniform along each analysesd set. As consequence of the above hypothesis, the standardized errors ε * o , ε * b , are also independent of σ cb . Using the above definitions, the innovations can be written as:

d = σ co ε * o -σ cb ε * b + b co -b cb , (7)
in which d is a function of two random variables: ε * o , ε * b , and the background-error standard deviation σ cb , running through all the possible tuned values σ cb , produced by the data assimilation system. Under the above hypotheses, the average, the centred moments and the cumulants of d are then factorized in terms of assumed single and joint standardized error moments, as well as in terms of power averages of σ cb .

Therefore,the innovation bias then comes as:

d ≡ b d = σ co E (ε * o ) -σ cb E (ε * b ) + b co -b cb = b co -b cb , ( 8
)
where E is the expectation operator taken over the standardized errors and the bar operator runs through the sample. ¿From that expression, we get the departure d = σ co ε * o -σ cb ε * b , from which centred moments of the innovation are computed. Taking into account that

E(ε * 2 b ) = E(ε * 2 o ) = 1
, by definition, one writes the variance of innovation as:

d 2 ≡ σ 2 d = σ 2 co + σ 2 cb -2σ co σ cb c ob , (9)
where

c ob = E(ε * o ε * b )
is the error correlation. Under the usual hypothesis of uncorrelated errors, the innovation variance is simply the sum of the mean error variances:

σ 2 o ≡ σ 2 co ≡ β o σ 2 d ; (10a) σ 2 b ≡ σ 2 cb ≡ β b σ 2 d . (10b)
where β o and β b are the explained variance fractions of the innovation variance. Uncorrelated errors lead to:

β o + β b = 1.

Skewness and Kurtosis

The formulae of the skewness s d and kurtosis k d of innovations are obtained after expanding the third and fourth order powers of d and then averaging it. For that purpose we need to define the single error skewnesses:

s co ≡ E ε * 3 o ; ( 11a
)
s cb ≡ E ε * 3 b , (11b)
the single error kurtosis:

k co ≡ E ε * 4 o -3; (12a) k cb ≡ E ε * 4 b -3, (12b)
the joint cumulants of order 3:

κ oob ≡ E (ε * o ε * o ε * b ) ; (13a) κ bbo ≡ E (ε * b ε * b ε * o ) (13b)
and the joint standardized error cumulants of order 4:

κ ooob ≡ E (ε * o ε * o ε * o ε * b ) -3E (ε * o ε * b ) ; (14a) κ oobb ≡ E (ε * o ε * o ε * b ε * b ) -2E (ε * o ε * b ) 2 -1; (14b) κ bbbo ≡ E (ε * b ε * b ε * b ε * o ) -3E (ε * o ε * b ) . (14c)
If observation and background errors are assumed bivariate Gaussian, then all the above, single and joint cumulants will vanish. If errors are independent, then only the joint cumulants vanish. The joint cumulants appear also in correlations between nonlinear functions of errors. The writing of s d and k d depends also on the average of powers of the tuned background error standard deviation σ cb . We express those averages in a non-dimensional form as:

h nb ≡ σ n cb σ 2 cb n/2 = σ n sb σ 2 sb n/2 , n = 1, 2, 3, 4. (15)
The h-terms h nb work as measures of the heteroscedasticity of background errors. Since the tuning factor f b is assumed uniform, those terms are simply obtained from statistics of σ sb . The h-terms are equal to 1 if the specified σ sb does not vary along the sample, i.e. in homoscedastic conditions. By definition, h 2b = 1.

We still have the inequalities: h 1b ≤ 1, h 3b ≥ 1 and h 4b ≥ 1, thanks to the Jensen inequality [START_REF] Rudin | Real and Complex Analysis[END_REF], applied to the convex power functions of the type: The skewness s d of innovations is then written as:

g(u) = u p (u = σ 2 cb , p ≥ 1).
s d = β 3/2 o s co -β 3/2 b h 3b s cb + 3β b β 1/2 o κ bbo -3β o β 1/2 b h 1b κ oob , (16)
in terms of the skewnesses of errors s co , s cb , and of the third-order joint error cumulants, weighted by powers of the variance fractions β o and β b and heteroscedasticity h-terms. We verify that observation and background errors skew the innovations, towards respectively the same side and the opposite side.

Innovations can be skewed even when errors are uncorrelated (c ob = 0) and have Gaussian marginal pdfs (s co = s cb = 0). That is due to non-linear error correlations brought by third-order joint error cumulants.

Dividing Equation ( 9) by σ 2 d , one gets:

β o + β b -2c ob h 1b β 1/2 o β 1/2 b = 1. ( 17
)
Playing with the above equation, the kurtosis k d of the innovation comes as:

k d = k d1 + k d2 + k d3 ; (18a) k d1 = β 2 o k co + β 2 b h 4b k cb ; (18b)
k d2 = 6β o β b κ oobb - 4β 3/2 o β 1/2 b h 1b κ ooob -4β 1/2 o β 3/2 b h 3b κ bbbo ; (18c)
k d3 = 3β 2 b (h 4b -1) + 12c ob β 1/2 o β 1/2 b β b (h 1b -h 3b ) -β 1/2 o β 1/2 b c ob (h 1b -1) 2 . ( 18d
)
The term within k d1 in Equation (18a) is the contribution from the kurtosis of individual errors. The second term k d2 is uniquely due to the fourth-order standardized error joint cumulants and finally the third term k d3 is due to heteroscedasticity. The kurtosis k d can be completely attributed to one of the above terms. This means that a non-vanishing k d can occur for non-independent, uncorrelated errors with Gaussian marginal pdfs (k co = k cb = c ob = 0). On the other hand, bivariate Gaussian independent errors in a sample of heteroscedastic background errors produce a symmetric leptokurtic innovation pdf with s d = 0

and

k d = 3β 2 b (h 4b -1) > 0.
This is well understood because the mixture of Gaussian random variables of different variances is non-Gaussian.

Constraints under independent errors

The statistical error independency is a common hypothesis taken in data assimilation. Let us present here some of the constraints of error skewnesses and kurtosis imposed by error independency. The averaged skewness s b and averaged kurtosis k b of the background error over a heteroscedastic sample are respectively given by:

s b = h 3b s cb ; (19a)
k b = h 4b k cb + 3 (h 4b -1) . (19b)
Since observational errors are homoscedastic, we get directly the mean skewness and the mean kurtosis:

s o = s co ; (20a) k o = k co . (20b)
When errors are assumed statistically independent, the skewness s d and the kurtosis k d are given as linear combinations, respectively of the error skewnesses: (s o , s b ) and of the error kurtosis (k o , k b ):

s d = β 3/2 o s o -β 3/2 b s b ; (21a)
k d = β 2 o k o + β 2 b k b . (21b)
The corresponding weights are proportional to powers of error variance fractions β o and β b suggesting that non-Gaussianity of innovations tends to be explained by the non-Gaussian error with the largest variance.

Under independent errors, the value of η d in (3) has a stricter lower bound provided by the positiveness of the following variance:

var r d + α ε o ε b σ 2 d = η d + β o β b α (α -4) ≥ η d -4β o β b ≥ 0, ∀α, (22)
with r d meaning the same as in Equation ( 3). This sets a necessary condition of error's independency for fixed variance fractions of errors. In all the 13 analysed samples, η ). This is illustrated in the example of Fig. 3 by delimiting the possible domains of (s b , k b ) (Fig. 3a) and (s o , k o ) (Fig. 3b), compatible with conditions:

s 2 b -2 ≤ k b ≤ β -2 b k d + 2β 2 o -β b β -1 o s b + β -3/2 b s d 2 ; ( 23a
)
s 2 o -2 ≤ k o ≤ β -2 o k d + 2β 2 b -β o β -1 b s o -β -3/2 o s d 2 . ( 23b
β o = 0.25,β b = 0.75, s d = 1 and k d = 1.
In particular, under the above conditions, the background error, the one with largest variance, is not allowed to be Gaussian, since the pair (s b = 0, k b = 0) is excluded from the allowed set as seen in Fig. 3a. This suggests that errors with largest variances are forbidden to be Gaussian as long as consistence with non-Gaussianity of innovations is verified. When heteroscedasticity is present, the skewnesses and kurtoses must be substituted by (19a), (19b), (20a), (20b) in the inequalities (23a), (23b).

We have analyzed the conditions for the full attribution of non-Gaussianity to one of the errors, keeping the remaining one to be Gaussian. Playing with inequalities (23a), (23b), and taking (s o = 0, k o = 0) and/or (s b = 0, k b = 0), we prove that the necessary condition for the existence of a Gaussian error of variance fraction β (i.e β o and/or β b ) is given by the following third-order polynomial inequality:

P (β) ≡ β 3 -3β 2 + β 3 + k d 2 - η d 2 ≤ 0 ; β ∈ [0, 1] . (24)
The variance fraction of the non-Gaussian error is 1 -β. The polynomial P (β) has a single root when the limiting curve η d = 0 is approached thus requiring increasing variance to the non-Gaussian error with the variance fraction 1 -β g . In practice, given specified consistent values β o and β b and comparing them to β g , it is possible to state whether non-Gaussianity of innovations can exclusively come from one of the errors (scenario of Single Source of non-Gaussianity) in alternative to its sharing among errors (scenario of Multiple Source of non-Gaussianity). These two scenarios will be explored in sections 5 and 6 for the analysed samples. The reverse condition of 24 states that if β is the variance fraction of a certain error and P (β) > 0, then that error cannot be Gaussian or the errors are not independent.

β g in the interval [0, 1] which is function of (s d , k d ). Since P (0) = -η d /2 ≤ 0 and P (1) = s 2 d ≥ 0,
The consistent tuned error fractions β o , β b depend on tuning factors f o , f b . This means that the partition of innovation's non-Gaussianity among errors is not independent from the problem of variance error tuning.

Therefore bounds are imposed to tuning factors in the case of single attribution of non-Gaussianity. Let β so (β sb ) be the quotient between the mean of specified observation (background) error variances and the innovation variance. Then the necessary conditions for Gaussian observation errors or equivalently the necessary condition for non-Gaussianity coming exclusively from background errors is:

f o ≤ (β g /β so ) 1/2 ≡ f o1 ; (25a) β -1/2 sb ≥ f b ≥ [(1 -β g ) /β sb ] 1/2 ≡ f b2 (25b)
On the other hand, the necessary conditions for Gaussian background errors or equivalently the necessary conditions for non-Gaussianity coming exclusively from observation errors is:

f b ≤ (β g /β sb ) 1/2 ≡ f b1 ; (26a) β -1/2 so ≥ f o ≥ [(1 -β g ) /β so ] 1/2 ≡ f o2 (26b)
The values of β g and the limiting tuning factors f o1 , f o2, , f b1 , f b2 are listed in Table 1 for the 13 studied samples. For example, under the hypothesis of Gaussian observation errors, that requires, for most of the samples, a maximum tuning f o1 of the order of 0.5, thus strongly deflating the specified variance of the observation error. We also notice the necessary inflation of background error variance for sample 4s at which f b2 = 1.305. Under the hypothesis of Gaussian background errors, we notice in particular the need for a strong deflation of the specified background error variances in the 5 'ice-covered' samples, as it is corroborated by the values f b1 lying in the interval [0.16, 0.54].

Error pdfs obtained by the Maximum Entropy Method

General Formulation

The application of the Maximum Entropy (ME) method is quite useful in data assimilation and inversion problems when pdfs of errors are non-Gaussian. Some applications of the ME in this context can be found in [START_REF] Urban | Humidity profile retrieval using a maximum entropy principle[END_REF][START_REF] Bocquet | Reconstruction of an atmospheric tracer source using the principle of maximum entropy I[END_REF][START_REF] Bocquet | Reconstruction of an atmospheric tracer source using the principle of maximum entropy II[END_REF][START_REF] Bocquet | Inverse modelling of atmospheric tracers: non-gaussian methods and second-order sensitivity analysis[END_REF]. In this section we write down explicit pdfs of errors that are consistent both with: a) the innovation statistical moments of orders one up to four and: b) the assumed marginal moments of errors.

Let us consider the joint pdf ρ ob of errors in the observation space and ρ o * b * the corresponding joint pdf of the standardized errors introduced by Equations (6a), (6b). Those pdfs are related by:

ρ ob (ε o , ε b ) = ρ o * b * (ε * o , ε * b ) σ co σ cb . ( 27
)
The least committing method of finding ρ ob (ε o , ε b ) is the Maximum Entropy (ME) method [START_REF] Jaynes | Information statistics and statistical mechanics i[END_REF][START_REF] Jaynes | Information statistics and statistical mechanics ii[END_REF][START_REF] Cover | Elements of Information Theory[END_REF],

looking for the maximum of the differential joint Shannon entropy [START_REF] Shannon | The mathematical theory of communication[END_REF] of errors H(ε o , ε b ), subject to a certain set of constraints. Each set of innovations, presented before, correspond to samples of realizations of the error pairs (ε o , ε b ). Then the pdf of errors is estimated by maximizing the entropy over the each set. Assuming that all those pairs are independent from each other due to the data thinning performed previously and also that standardized errors have the same pdf over the sample, we have the sampling mean Shannon entropy written as:

H (ε o , ε b ) ≡ E [-ln ρ ob (ε o , ε b )] = ln (σ co ) + ln (σ cb ) + H (ε * o , ε * b ) . ( 28
)
where

H(ε * o , ε * b )
is the Shannon entropy of the standardized errors. The problem of entropy maximization reduces to the search of

H M (ε * o , ε * b ) ≡ max[H(ε * o , ε * b )],
over the pdfs satisfying a set Ω of N c strong constraints in terms of expectations of functions f 1 ,...,f N c of the standardized errors:

Ω = {E [f k (ε * o , ε * b )] = c k ; k = 1, ..., N c } . ( 29
)
When H M is finite, it is found at the global minimum of a function Γ, globally convex in terms of the pdf, of N c Lagrange multipliers (LMs) λ 1 , . . . , λ N c , spanning the minimization control space:

H M (ε * o , ε * b ) = min λ1,...,λ N c , Γ (λ 1 , . . . , λ N c ) ≡ Γ λ1 , . . . , λNc , (30) with:
Γ (λ 1 ...λ Nc ) ≡ ln Z - Nc k=1 λ k c k ; Z (λ 1 . . . λ Nc ) ≡ exp Nc k=1 (λ k f k (ε * o , ε * b )) dε * o dε * b , ( 31
)
where Z is the partition or normalization function. The pdf of maximum entropy (ME-pdf) belongs to the exponential family [START_REF] Barndorff-Nielsen | Information and exponential families in statistical theory[END_REF] and reads:

ρ o * b * (ε * o , ε * b ) ≡ exp Nc k=1 λk f k (ε * o , ε * b ) Z λ1 , . . . λNc , (32)
where we have used the Lagrange multipliers found at the maximum entropy solution or minimum of Γ.

In general, fewer is the number N c of constraints, larger is the value of H M . When the set Ω includes only error moments of order one and two, then the ME-pdf is Gaussian. In our applications below, the

f k (k = 1, . . . , N c
) functions are polynomials of the errors, corresponding to error and innovation moments.

Maximum-Entropy errors constrained by Innovation statistics

We will restrict our study to the case of statistically independent errors. In this case the joint entropy of standardized errors is the sum of individual entropies: 28. In our case, the set of constraint functions restrict to sums of individual observation and background error monomials, thus leading to separable single error integrals in the formulation of the log-partition function ln(Z). Two kinds of maximum-entropy problems are formulated bellow. The first one is the so-called problem of 'Prescribed error non-Gaussianity'. Here, each observation and background standardized errors satisfy the N c = 4 constraints:

H(ε * o ) + H(ε * b ) in Equation
E [ε * t ] = 0; E (ε * t ) 2 = 1; E (ε * t ) 3 = s ct ; E (ε * t ) 4 = k ct + 3; t = o, b, (33)
with s co , s cb , k co , k cb , being defined by Equations (11a), (11b), (12a), (12b), and satisfying innovation constraints given by Equations ( 16), (18a). A particular case is that in which all non-Gaussianity is fully attributed to one of the errors i.e. the case of 'Single source of non-Gaussianity'. The (s, k) or (skewness, kurtosis) dependence of the maximum entropy H M and of the four Lagrange multipliers (LM) λ p (p=1,..,4), corresponding to the powers p in Equation 33 is discussed by [START_REF] Wu | Calculation of maximum entropy densities with application to income distribution[END_REF] as well in the appendix.

The second kind of ME problem, hereafter called 'Maximum entropy sharing of non-Gaussianity', correspond to the case in which the skewnesses and kurtoses of errors are not anchored to prescribed values, though they still satisfy the innovation constraints given by Equations (21a), (21b). In this case we get a total of N c = 6 constraints terms in the standardized errors:

E [ε * t ] = 0; E (ε * t ) 2 = 1, with t = o, b; (34a)
E ε * 3 o β 3/2 o -ε * 3 b h 3b β 3/2 b = s d ; (34b) E ε * 4 o β 2 o + ε * 4 b h 4b β 2 b = k d + 3 β 2 o + β 2 b , (34c)
to which correspond 6 Lagrange multipliers: λ o1 , λ o2 , λ b1 , λ b2 , λ d1 , λ d2 with the first and second pair of LMs, corresponding to the mean and variance constraints, respectively for the observation and background errors. The last two LMs correspond to constraints related to s d and k d respectively. The function Γ to be minimized reads:

Γ (λ o1 , λ o2 , λ b1 , λ b2 , λ d3 , λ d4 ) = t=b,o ln exp λ t1 ε + λ t2 ε 2 dε + ln exp λ d3 β 3/2 o ε 3 + λ d4 β 2 o ε 4 dε + ln exp -λ d3 β 3/2 b h 3b ε 3 + λ d4 β 2 b h 4b ε 4 dε - λ b2 + λ o2 + s d λ d3 + λ d4 k d + 3 β 2 o + β 2 b . (35)
The corresponding ME-pdfs of the standardized errors are written in the condensed equation:

ρ t * (ε) = exp (Q(ε)) exp (Q(u)) du (36a) Q(ε) = λt1 ε + λt2 ε 2 + λd4 β 2 t h 4t ε 4 ± λd3 β 3/2 t h 3t ε 3 (36b)
Positive and negative signs in (36b) apply respectively to the observation error (t = o) and to the background error (t = b) while using h 3o = h 4o = 1 in the equation of the pdf of ε * o . The LMs are those minimizing the function [START_REF] Gilbert | Some numerical experiments with variable storage quasi-newton algorithms[END_REF].

For background homoscedastic errors, we have: by minimizing the approximate error statistical negentropies, defined as in Equation ( 5) and subject to constraints given by Equations (21a), (21b), thus leading to:

h 3b = h 4b = 1
(s o , s b ) ≈ s d β 3 o + β 3 b β 3/2 o , -β 3/2 b ; (37a) (k o , k b ) ≈ k d β 4 o + β 4 b β 2 o , β 2 b . (37b)
As a consequence, when one of the errors is dominant, i.e. max(β o , β b ) is close to unity, the non-Gaussianity of innovations through (s d , k d ) tends to be explained by the error with the largest variance, whereas the remaining one is close to Gaussianity.

We will refer again to the example of Figs. All of the above ME-problems are numerically solved following a modified version of the algorithm presented in [START_REF] Rockinger | Entropy densities with an application to autoregressive conditional skewness and kurtosis[END_REF], also used in [START_REF] Pires | Non-gaussianity and asymmetry of the winter monthly precipitation estimation from nao[END_REF]. The support of ME-pdfs is the finite interval [-20, 20], which seems satisfactory in our applications. Numerical details are given in the appendix.

Minimum Variance and Best Linear Unbiased Estimators

Generic Properties

When observation and background errors are statistically independent, as previously assumed, then the Minimum Variance Unbiased Estimator (MVUE) x M , issued from an observation y and a background estimate x b = y -d is given through the convolution of the marginal error pdfs in the following form:

x M = x b -b cb + ∆ M b = x + ε M ; ε M = ε b + ∆ M b ; (38a) ∆ M b ≡ -E (ε b |d) = -σ cb uρ b * (u) ρ o * (d -b d + uσ cb ) σ -1 co du ρ b * (u) ρ o * (d -b d + uσ cb ) σ -1 co du , (38b)
where the ME-pdfs of errors habe been used and x is the true state value. An alternate estimator to the MVUE is the Maximum A Posteriori (MAP) corresponding to the absolute mode of the a posteriori pdf.

The error of the MVUE is ε M = x -x M and ∆ M b is the correction to be added to the unbiased background in order to obtain MVUE. We notice that ∆ M b is the symmetric of the conditional mean of ε b for a given innovation value d. The normalization integral in denominator of Equation (38b) is proportional to the estimated conditional-innovation theoretical pdf ρ d (d|σ cb ), inferred from ME-pdfs of errors and computed at the particular value of d and σ cb . Then the innovation pdf over a certain set must be approached by the mean of ρ d (d|σ cb ) running over all the innovations of that set. This pdf is comparable to the the empirical raw innovation pdf, thus providing another consistency test beyond those given by the skewness and kurtosis, as discussed in Section 3. However, the mean theoretical pdf does not necessarily coincide with the ME-pdf constrained by innovation statistics only.

The mean square of the MVUE error is given by the conditional variance of ε b computed at the innovation value d, i.e. E(ε 2 M |d) = var(ε b |d). For particular extreme innovations, say larger than two standard deviations and certain non-Gaussian error pdfs (e.g. exponential pdf laws), that posterior variance can be larger than the smallest of the error a priori variances i.e. the min(σ 2 b , σ 2 o ), [START_REF] Talagrand | Bayesian estimation. Optimal interpolation. Statistical linear estimation[END_REF]. The same inequality holds for the conditional a posteriori Shannon entropy [START_REF] Pires | Aplicação da entropia estatística à assimilação de dados[END_REF]. In spite of that, the overall mean-square of the a posteriori error, say for a large sample of innovations values d, shall be smaller than a priori mean square error: min(σ 2 b , σ 2 o ). In fact, this is due to the Law of Total Variance in the form:

E d ε 2 M ≡ σ 2 M = E d E ε 2 M |d = σ 2 b -E d E (ε b |d) 2 ≤ min(σ 2 b , σ 2 o ), (39)
where the expectation E d is the average running through a large sample of innovations. When errors are Gaussian, the MVUE is the classical Best Linear Unbiased Estimator (BLUE), [START_REF] Gelb | Applied Optimal Estimation[END_REF] given by:

x L = x b -b cb + ∆ Lb = x + ε L ; ε L = ε L + ∆ Lb ; (40a) ∆ Lb ≡ σ 2 cb σ 2 cb + σ 2 co (y -x b -b d ) , (40b)
where ∆ Lb is the linear correction to be added to the unbiased background and ε L is the BLUE error.

Even when errors are non-Gaussian, the BLUE still gives the optimal linear estimator of the state x from the observation y and the background x b . The impact in the a posteriori estimate of taking into account non-Gaussianity of errors is given by the difference between the BLUE and MVUE:

∆ M L ≡ x M -x L = ∆ M b -∆ Lb = ε M -ε L , (41)
also given by the difference between the BLUE and MVUE errors. The mean square error of the BLUE, over a large set of innovations, admits the decomposition:

E d ε 2 L ≡ σ 2 L = σ 2 d β o β b = σ 2 M + E d ∆ 2 M L ≥ σ 2 M , (42)
which shows that a self-consistent assimilation scheme based on the non-Gaussian MVUE must in general perform better than that based on the BLUE.

Average skill in the SK-domain

In this section we assess the performance of the MVUE over the innovation SK-domain. In order to evaluate the overall relative skill of the MVUE, compared to the BLUE, we compute the non-dimensional Mean Square Error Skill Score (MSESS), taken with respect to the BLUE, also often used in forecasting evaluation [START_REF] Déqué | Continuous variables, chapter 5[END_REF]:

MSESS ≡ E d ∆ 2 M L σ 2 L = 1 - σ 2 M σ 2 L ∈ [0, 1] . (43)
This is independent from the innovation variance and ranges between 0, for Gaussian errors and 1 when the MVUE is exact, i.e. null a posteriori error, holding if one of the error variances vanishes. The MSESS depends on the error variance partition: β o , β b = 1 -β o and also on the assumed pdfs of the standardized errors. The MSESS can also be seen as a measure of the BLUE sub-optimality in presence of non-Gaussian errors.

Next, we study how MSESS changes with the fixed constraints for the three scenarios of non-Gaussianity partition among errors, described in section 5.2. Finally, values of the MSESS are assessed for the analysed samples of brightness temperatures.

Single-source of non-Gaussianity

We assess the MSESS, defined by Equation (43) when only one of the errors is non-Gaussian. Without loss of generality we will attribute all non-Gaussianity to the background error of variance β b for which:

(s b , k b ) = (-s d /β

Prescribed non-Gaussian errors

In a second experiment, we study the case of two non-Gaussian errors. In order to mimic that situation, we compute MSESS by fixing equal error variances: β o = β b = 0.5). For the observation error we fix

Maximum entropy sharing of non-Gaussianity

Here, we assess how the MSESS depends on both innovation non-Gaussianity through (s d , k d ) and the error variance partition (β o , β b ) under the scenario of Maximum entropy sharing of non-Gaussianity.

From symmetry considerations, the MSESS is symmetric with respect to s d and invariant for the mapping

β o ↔ 1 -β o ,

Average skill for brightness temperature innovation samples

One of the aims of this article is to assess the potential average gains from using the MVUE, inferred from the Maximum Entropy method, as an alternative to the BLUE. This is applied here to the assimilation of brightness temperatures extracted from HIRS channels. In order to measure that, we use an adimensional score (larger when it is better), based on the mean square of the correction (MVUE minus BLUE): ∆ M L =

x M -x L , computed for different classes of standardized innovation:

d * = (d -b d )/σ d over the innovation
samples. The score is obtained from conditional averages over three classes of the standardized innovations according to its size as:

MSC j = σ -2 d ∆ 2 M L | |d * | ∈ I j ; I a = [0, 1[; I b = [1, 2[; I c =]2, ∞[. ( 44
)
where MSC j (j = a, b, c) stands for Mean Score over Classes, where MSC a , MSC b , and MSC c , are computed for successively higher values of |d * |, within one, two and over two standard deviations. In the present study the analyzed innovation samples are not very far from Gaussianity. In fact, the pairs (s d , k d ) from Table 1 are quite close to point (0, 0) and, therefore, from

p j MSC j ; (45a) MSESS = MSC β o β b . ( 45b
)
The three scores MSC j (j = a, b, c) are computed in a quite exhaustive way for 12 brightness temperature innovation samples, where non-Gaussianity is statistically significant. The sample 14s has been excluded due to its non-significant non-Gaussianity. For each sample, we have considered the three following scenarios: Finally, we notice that only certain tuning intervals are compatible with the hypothesis of Single-source of non-Gaussianity. For example, in sets 4s, 5s, 7s and all 'ice-covered' samples, the absence of background error tuning, i.e. f b = 1 do not match with an exclusive non-Gaussian background error. This shows that tuning of error variances and the attribution of non-Gaussianity to errors must be integrated in a more general problem.

In order to see how the MVUE and the BLUE correct the background estimator, we have considered a case for the most skewed sample of innovations (4s), where s d = -0.70, k d = 1.02 and the heteroscedasticity is quite small (Cf. Fig. 2). The error variance partition is:

β b = 0.47, β o = 0.

Discussion and Conclusions

Histograms of innovations for certain assimilated observables (e.g. brightness temperature obtained from HIRS channels) reveal statistically significant non-Gaussianity, as far as skewness and kurtosis are concerned. Non-Gaussianity of innovations can be due to: marginal non-Gaussian observation or background errors, non-linear error correlations and ultimately to heteroscedasticity of errors within the sample used to estimate innovations statistics. The origin of non-Gaussianity can only be determined by independently tested hypothesis that would suggest the contributions of the instruments, observation representativeness, models and assimilation algorithms. When non-Gaussianity of innovations essentially comes from one of the errors, then its variance must be sufficiently high in order to accommodate it. This imposes some bounds on the error variance tuning. In some cases, it leads to incompatibility between absence of tuning and specific sources of non-Gaussianity, i.e. variance tuning is compulsory while taking into account non-Gaussianity in the data assimilation schemes.

The least committing method of taking into account non-Gaussianity in the error's pdfs is the Maximum Entropy method. The resulting pdfs are then used to get the Minimum Variance Unbiased Estimator (MVUE) of the system state in the observational space. We analyse which conditions favour the largest increment of performance when compared to the Best Linear Unbiased Estimator (BLUE), vastly used in data assimilation. This relative performance is measured here by the MSESS (Mean Square Error Skill Score), ranging between 0 (no added value) up to 1 (total correction). The MSESS is particularly large, i.e., the role of non-Gaussianity is high, when errors from a priori estimates (observations, backgrounds etc.), are skewed and in majority towards the same side. Compensation of error skewnesses degrades the performance of non-Gaussian data assimilation schemes. Moreover, in what concerns data errors and their non-Gaussianity impact: 'Large skewness put into a small variance' has a larger impact than 'Small skewness put into a large variance'. In the present paper, we have applied the Maximum Entropy method to compute error pdfs that are consistent with innovation statistics of the brightness temperature, obtained from HIRS channels. The largest difference between the MVUE estimator and the BLUE estimator holds at largest absolute innovations, typically larger than two standard deviations σ d of innovations and reaching 10% -60% of σ d . This can eventually be applied in operational data assimilation schemes. The first step would be the computation of analysis obtained by the MVUE in the observation space which in turn would be re-assimilated as corrected unbiased data and also with the effect of non-Gaussianity of errors taken into account.

The implementation of non-Gaussian error models into data assimilation is still in its birth. Many aspects still require new ideas such as how to generalize non-Gaussian error models to the whole innovation vector.

There, models for multivariate skewnesses and kurtoses [START_REF] Takafumi | On a measure of multivariate skewness and a test for multivariate normality[END_REF] have to be taken into account and also the effect of spatial error correlations. Another goal is how to state error's non-Gaussianity in the spectral domain, often used by numerical forecasting models and how to compute it from the ensemble methods along with the background error covariance matrices and cumulant tensors of order three and four. Another possible future extension is the addition of penalty terms to data assimilation cost-functions, weighting the deviation of error's pdfs from Gaussianity. A systematic similar study to innovations of other observable types is in mind.

Acknowledgements

This research was developed at IDL and LMD with support from the project: 'Développements méthodologiques

dH M (ε) dL = ρ (L) + ρ (-L) , ( 46
)
where ρ is the ME-pdf. The logarithm ln ρ at the interval bounds is scaled as O(-L 4 ), decreasing as far as L increases. That exponent, 4, is here the maximum even power of the constraining moments. Due to the finite L and for pdfs with kurtosis k and skewness s, near the limit k ≈ s 2 -2, a small hump in the pdf tail will accommodate the assumed innovation constraints. Given the scaling above and the range of moments, we can extend the pdf support D sufficiently further so as to get negligible bound effects, both on the entropy and on the ME-pdfs. Furthermore, in order to get integrands of order exp(O(1)), during the optimization process, we solve the ME problem for the scaled variable ε/L in the interval [-1, 1] by taking the appropriate scaled constraints. Afterwards, we apply the scaling entropy relationship:

H M (ε) = H M (ε/L) + log(L), (47)
where H M (ε/L) is the entropy of the scaled error. The LMs also depend on the interval amplitude L, but stabilizing for L high enough. After differentiating the expectation constraints for fixed moments c k of order k = 1, 2, 3, 4 = N c , we get:

dλ k dL = Nc=4 i=1 M -1 k,i × ρ (L) L i -c i + ρ (-L) (-L) i -c i ; k = 1, ..., 4, (48)
where M is a square matrix of order N c whose elements M k,i are the covariances between the powers ε k and ε i as determined by the ME-pdf for L. The elements M -1 k,i in (48) are the (k, i) elements of the inverse matrix of M. The minimization problem is solved by the Quasi-Newton method using the M1QN3 routine [START_REF] Gilbert | Some numerical experiments with variable storage quasi-newton algorithms[END_REF], starting at λ k = 0 for all k. The integrals giving the functional Γ to minimize and its λ-derivatives are approximated by the Gauss quadrature rule with N p weighting factors in the interval [-1,1]. In order to get full resolution during the minimization, and to avoid NAN (Not-a-Number) and INF (Infinite) numbers in computation, we subtract the polynomials in the arguments of exponentials by the correspondent maximum in D. Finally, the functional Γ is multiplied by a sufficient high factor F a in order to emphasize the gradient.

Considering the range of constraint moments, several experiments have led to the reasonable values: L = 20, N p = 500, F a = 1000, for which convergence is obtained nearly at 60 iterates with an accuracy of 10 -6 of the gradient of Γ. Larger values of L require larger values of N p .

Table 1: Statistics of the background errors, observation errors and innovations for the samples 4s, 5s, 6s, 7s, 11s, 12s, 14s 15s for 'ice-free' conditions and samples 5i, 6i, 7i, 14i and 15i for 'ice-covered' conditions. Table lists: See text for details.

  variable r d is the second order polynomial of the standardized innovation d * ≡ d /σ d , which has unitary quadratic coefficient and has the minimum variance. On the other hand, we easily get any pair (s d , k d ) satisfying (3) by using a pdf mixture of three weighted Dirac-Deltas. The set of pairs (skeweness, kurtosis) satisfying inequality (3) is hereafter called the SK-domain. The value of η d for a Gaussian pdf equals 2, whereas in the studied samples, η d ranges from 1.49 (sample 6i) to 3.99 (sample 12s).

≥ 1 ,

 1 and thus not imposing bounds to β o , β b since max(4β o β b ) = 1. From the Equations (21a), (21b) of s d and k d we can invert s o and k o , and/or s b and k b . Then applying inequalities of type (3) to the observation and background errors, i.e. k o -s 2 o + 2 ≥ 0 and k b -s 2 b + 2 ≥ 0 with the substituted s o , s b , k o and k b leads to the following bounds:

  ) Therefore, giving consistent error variance fractions β o , β b and innovation statistics (s d , k d ), there are bounded sets for the skewness and kurtosis of errors (k o , k b ) delimited by two parabolas of opposite concavity in the SK-domains of (s o , k o ) and (s b , k b ), i.e. those domains delimited by inequalities (23a), (23b). We easily verify that the necessary and sufficient condition for the non-vanishing of the set of solutions sets is given by the inequality (22). The biggest of the two sets of allowable pairs (s o , k o ) and (s b , k b ) corresponds to the error of smallest variance fraction (β o or β b

  one concludes that the maximum variance fraction of the Gaussian error is β g (s d , k d ) and the required minimum variance of the non-Gaussian error is 1 -β g (s d , k d ). In other words, (24) is equivalent to: β ≤ β g (s d , k d ). The contour field of β g (s d , k d ) appears in Fig. 4. The value of β g (s d , k d ) is symmetrical with respect to s d and decreases

  and (s b , k b ) = (s cb , k cb ). In this particular case, the ME solution is uniquely defined by the triplet (s d , k d , β o ). As consequence of the anti-symmetric roles of errors, the corresponding skewnesses s o and s b are anti-symmetrical with respect to s d and invariant in the mapping: (s o , β o ) ↔ (-s b , 1 -β o ). Furthermore, k o and k b are symmetrical with respect to s d and invariant for the mapping: (k o , β o ) ↔ (k b , 1 -β o ). The pairs (s o , s b ) and (k o , k b ) are approximately given

  3a, 3b, corresponding to (s d =1,k d =1) and β o =0.25, β b =0.75, in order to show where the ME-solution lies in the allowable SK domains of the observation and background errors. The contours of the sum of maximum individual entropies:H M = H M (ε * o ) + H M (ε * b ), are drawn as function, both of (s b , k b ) (Fig.3a) and of (s o , k o ) (Fig.3b). The field of H M has a minimum of ∼ 1.4 nats at the edge of the allowable SK domains. Nat is the natural logarithmic unit of Shannon entropy and information with the conversion 1 nat = log 2 e bit ≈ 1.44bit. The maximum of H M ∼ 2.55 nats, lies at the interior of the domain, at (s o , k o ) ∼ (0.65, 0.71) and (s b , k b ) ∼ (-1.41, 1.70) whereas Equations (37a),(37b), give the close estimation (s o , k o ) ∼ (0.29, 0.19) and (s b , k b ) ∼ (-1.48, 1.75).

3 / 2 b

 32 , k d /β 2 b ). Conversely, the observation error has variance β o = 1 -β b . The MSESS contours for the case of equal error variances: β o = β b = 0.5 are shown in Fig. 5a for |s b | < 4.5 and -2 < k b < 17. The MSESS increases monotonically with |s b |, for any k b , from the zero value, at the semi-straight line: s b = 0, k b ≥ 0 (dashed in Fig. 5a), corresponding to leptokurtic non-Gaussian errors, up to a maximum value near the limit curve: η b ≡ k b -s 2 b + 2 = 0. There, the non-Gaussian error pdf is tightly constrained, leading to a singular behaviour of the posterior pdf and to a concentration of MSESS contour lines. The MSESS grows with the absolute value |k b | for platykurtic non-Gaussian errors (k b < 0), contrary to the case of leptokurtic errors (k b > 0) where MSESS slightly decreases with kurtosis and fixed skewness. Similar qualitative MSESS behaviour holds for different error variance partitions. Still keeping the same scenario, we study how MSESS changes with the Gaussian error variance β o and the skewness s b by fixing the non-Gaussian background kurtosis. The corresponding MSESS contours are shown in Fig. 5b for the particular case k b = 3. There, the MSESS field is symmetrical with respect to s b , which occurs in general when one of the errors is un-skewed i.e. s o = 0. As expected, the MSESS tends to zero (exact estimation) when β o tends to 0 or 1, corresponding respectively to exact observation and exacto background. Then, for a given s b , the maximum MSESS lies at a particular value of β o , close to 0.5, corresponding to the error variance partition for which the a posteriori BLUE variance is maximum. The value of Gaussian error variance β o (the Non-Gaussian error variance β b ) at the MSESS maximum, slightly decreases (increases) as far as the error absolute skewness increases. When the largest error variance comes from the Gaussian error (high β o ), the BLUE is rather satisfactory as it is seen from the low MSESS values (Fig. 5b). On the other hand, if the dominant error is non-Gaussian (low β o ), the MVUE is much better than the BLUE, especially at high values of |s b |. Now, fixing k d and s d and increasing β o up to the allowable maximum β g (s d , k d ) leads to increasing values of the skewness and kurtosis of the non-Gaussian error and also to an increasing MSESS (Figs. 6a-c) for k d = 1, 5, 13 respectively). This means that the smaller is the variance of the error to which non-Gaussianity is attributed, the larger is the difference between the MVUE and BLUE. If the non-Gaussian error variance β b is much lower β o much closer to 1), then k b can assume extremely high values with remarked effects of the finiteness of the ME-pdfs supports (e.g. k d = 5 and β b = 0.1 leads to k b = 503). This effect happens inside hill-shaped regions bounded by contour inflections in Figs. 6a-c.

(

  Fig. 5b for a fixed k b = 3. There, MSESS increase when one approaches the limits s b = ±(k b + 2) 1/2 . For s b ,s o of the same sign (right side of Fig. 7b), the MSESS grows regularly up to ∼ 0.8, whereas in the left side (opposite signs) it hardly reaches 0.2 or at much 0.5 for high values of |s b |. When errors are skewed towards the same side, i.e. for fixed s b > 0 (right side of Fig. 7b), the MSESS is maximum for an error variance partition such that β o -β b is a monotonous increasing function of |s o -s b |. In other words, high MSESS is favoured by high (small) error's skewness weighted by high (small) variances. On the other hand, if s b ,s o have opposite signs (left side of Fig. 7b), high MSESS holds when the error of largest variance has the largest absolute skewness. Equal error variances under opposite error skewnesses tend to produce small differences between the MVUE and BLUE. Although errors have finite non-null variance, the effect of error non-Gaussianity can have effects ranging from nearly exactness of the MVUE i.e. MSESS=1 up to non-existence of differences between the MVUE estimator and the BLUE (MSESS=0). This is illustrated in MSESS contouring of Fig. 3c in the allowable SK-domain of observation errors for the case (s d = 1, k d = 1, β o = 0.25). The MSESS field ranges between the maximum ∼ 0.98, lying near the left edge of the permitted SK domain. The minimum of MSESS, ∼ 0.01 lies at (s o , k o ) ∼ (2.37, 5.94), close to the right edge of the set.

2 d - 2 .

 22 due to the symmetrical roles of errors. The ME-solution, maximizing entropy, is not necessarily the one minimizing MSESS as it is clear for the case discussed in Fig. 3c. In fact, at the ME solution H M ∼ 2.55 nats, (s o , k o ) ∼ (0.65, 0.71), and the MSESS is ∼ 0.14. At the point of a Gaussian error, H M ∼ 2.45 nats, (s o , k o ) = (0, 0), and the MSESS is ∼ 0.25. Nevertheless, the value of MSESS is not in general a monotonous decreasing function of H M . Contours of MSESS over the (s d , k d ) plane are shown in Figs. 8a, 8b, for β o = 0.25 and 0.4 respectively. We notice, in the above figures, that k d has a lower bound given by Equation (22) in which k d ≥ 4β o (1-β o )+ s By setting high enough values of β o and k d , the MSESS increases from nearly zero at s d = 0, where the skewnesses of errors compensate themselves, and then reaching a maximum at a specific intermediate |s d | value and then decreases towards the SK domain edge (Figs. 8a,8b ). In Figs. 8c-e, we plot the MSESS as a function of |s d | and β o , respectively for k d = 1, 5, 13. We see that in the case of leptokurtic innovations (k d > 0) and relatively low skewness values, say |s d | 1, the MVUE and BLUE do not differ much, according to the values MSESS 0.01. The highest MSESS values hold when one of the error variances β o or β b ) range between 0.1 and 0.3, and the innovation skewness is sufficiently high. The Figures 9a-e show how the s o and k o change for different values of s d (in abscissas) and β o (in ordinates) for the same previous values of k d . Figs. 9a,c,e refer to s o and Figs. 9b,d,f refer to k o respectively for k d = 1, 5, 13. From those figures, one verifies that the ME solution tends to concentrate the skewness and kurtosis of the innovations onto the error with largest variance.

  Fig. 5a a small MSESS is expected as compared to the cases of extreme non-Gaussianity shown in previous sections. Therefore, the relative frequencies p a , p b and p c of the classes I a , I b and I c defined, in (44), must be close to those given by the standard Gaussian, i.e. p a ∼ 68%, p b ∼ 26% and p c ∼ 6%, corresponding to conditional averages of |d * |: q a ∼ 0.4, q b ∼ 1.5 and q c ∼ 2.5 respectively for classes I a , I b and I c . The MSC j (j = a, b, c) is measured in terms of the innovation variance. Through (43), one shows that the average MSC and the MSESS values are given as: MSC = j=a,b,c

1 )

 1 Figures Aa-l, as functions of f b , or equivalently as functions of the variance fraction β b of the background error. The scores are computed for the above referred scenarios: 1, 2 and 3 and the three classes a, b and c each. The tuning values range within estimated bounds given by (25),[START_REF] Shannon | The mathematical theory of communication[END_REF] for the scenarios 1 and 2. All in all, there are 9 score graphics per sample in each of Figure Aa-l. The root mean square of the differences between the BLUE and MVUE are obtained by taking the root square of MSC j (j = a, b, c) and then multiplying it by σ d from Table1. Some conclusions, common to every one of the 12 samples can be drawn. Firstly, for fixed error variances (β o , β b ) and all the three scenarios of non-Gaussianity (1, 2 or 3), one verifies that the partial scores MSC j (j = a, b, c) tend be higher for the most extreme values of the absolute innovation |d|, i.e. MSC a < MSC b < MSC c . For the scenarios of Single Source of non-Gaussianity (scenarios 1 and 2) and certain samples, the partial values of MSC reaches values up to 0.1-0.6 for the extreme class I c , i.e. when |d| differs more than two standard deviations, corresponding to ∼ 6% of the cases. This is particularly evident in some 'ice-free' samples namely: 4s, 5s, 6s, 11s, and 15s in which MSC c , corresponding to the most extreme innovation values, go beyond the threshold of 10%, or a mean correction for non-Gaussianity ∆ M L (MVUE minus BLUE) of about ∼ 0.31σ d , for a quite large range of possible background error variances, thus showing that, independently of the correct tunings, the relative performance of the MVUE is meaningful. Furthermore, in most of those cases, the variance of the non-Gaussian errors does not need to reach extremely low variances corresponding to unlikely high values of skewness and kurtosis. For the 'ice-covered' samples, the impact of non-Gaussianity in the MVUE is not very high.Secondly, in the medium classes I a and I b , the correction ∆ M L is quite small for the majority of consistent combinations of error skewnesses and kurtoses. For each scenario and class, the behaviour of MSC j (j = a, b, c) as a function of the variance of errors follows closely the behaviour described in section 6.2 in terms of MSESS. In fact, in scenarios 1 and 2 (Single-source of non-Gaussianity), MSC increases as far as the variance of the non-Gaussian error decreases, as suggested from Figs. 6a to 6c. There are quite local exceptions in samples 7s and 12s. In the scenario 3 (Maximum entropy sharing of non-Gaussianity), the partial values of the MSC present the lowest values, except in samples 4s and 5s where it reaches values close to 0.1. In that scenario, the partial scores reach local maxima at intermediate values of the error variance fraction, one 21

  Figure 11(b), we plot the corrections added to the background estimate due to the BLUE (∆ Lb ) and due to the MVUE (∆ M b ). The BLUE correction is linear in terms of the innovation d, with a weight given according to Equation (40b). The spread of points in both quantities is due to varying σ sb aling the set. From the definition of innovation d = ε o -ε b , positive (negative) values of d are favoured by negative (positive) values of ε b , which in our case are quite reduced (enhanced) due to the positive skewness of ε b . Therefore, as expected, for negative extreme values of the innovation the MVUE-correction: ∆ M b is quite enhanced, whereas for positive d a reduced correction is expected as shown in Fig. 11(b). In general, under the scenario of Single-source of non-Gaussianity, the MVUE-correction to the background tend to be enhanced with respect to the BLUE correction for extreme innovations d such that sgn(d) = sgn(σ d ). For a certain value d, a relatively marginal correction holds when: sgn(d) = -sgn(σ d ).

  pour l'assimilation des observations: Validation a posteriori, diagnostics d'optimalité et de Gaussianité' from the Programme National: 'Les Enveloppes Fluides et l'Environnement'. Thanks are due to Elias Holm for providing ECMWF feedback data assimilation files, Pedro Viterbo and Rui Perdigão for their comments and help. A. Appendix -Maximum entropy pdfs constrained by skewness and kurtosis In order to study the change of the Lagrange multipliers (LMs) and their ME-pdfs, we have solved the ME problem in the SK-domain for skewness s ∈ [0, 4.5] and kurtosis k ∈ [-2, 17] with steps of 0.05 in each parameter, using the pdf support set D = [-L, L] = [-20, 20] for numerical computations. In Figs.12a-ewe plot the contour lines, respectively of the Lagrange multipliers (LMs): λ 1 , . . . , λ 4 and of the Shannon entropy, as function of the imposed skewness and kurtosis. The LMs correspond to successive moments (p = 1, 2, 3, 4). In Fig.12f, we show the ME-pdfs for some typical pairs (s, k). The λ-fields are smoothly varying over the SK-domain except when they approach the boundary curve k = s 2 -2 where both the absolute values of the LMs and the negentropy tend to increase. There, the ME-pdfs exhibit very sharp picks. For symmetry reasons, the odd order λ-fields are anti-symmetric with respect to the axis s = 0 and the even order ones are symmetric. There are singular (s, k) points in which the contour lines of certain λ-fields tend to merge, as for example (s, k) = (0, 1) for λ 1 and λ 3 and (s, k) ∼ (1.5, 0) for λ 2 . The ME field, in Fig.12e, gets the maximum H G = 0.5 ln(2πe) nats at the Gaussian point (s, k) = (0, 0), decreasing from it as far the parabolic limit curve k = s 2 -2 is approached, both in the skewness and the kurtosis axes. When the support half-extent L → ∞, there are no ME-pdfs with s = 0 and k > 0[START_REF] Ormoneit | An efficient algorithm to compute maximum entropy densities[END_REF]. As consequence of finite L, λ 2 ∼ -0.5 and λ 1 ∼ λ 3 ∼ λ 4 ∼ 0 at the semi-straight line (s = 0, k > 0) as seen in Figs.12a-e. In spite of that, there are leptokurtic unskewed exponential pdfs (s = 0, k > 0) of the type ρ(ε) = C exp(-|ε| α ), α ≥ 2, not satisfying the ME requirement. In fact, at s = 0, the negentropy approximation of a ME-pdfs, given in Equation (5), is only a good approximation for k < 0.We have plotted in Fig.12f, the ME-pdfs at four typical points in the SK domain: G(Gaussian): (s, k) = (0, 0); SM(Skewed-Mesokurtic): (s, k) = (1, 1); SP(Skewed-Platykurtic): (s, k) = (1, -1); SL(Skewed Leptokurtic): (s, k) =(3.5, 12). The pdf of SM is very close to Gaussian as well as those with small s and high kurtosis k, since H M is very close to H G (Fig.12e). The pdf of SP is platykurtic with two modes and lower density near the origin. The pdf of SL is spiky near the origin with an extended tail of positive extreme anomalies and a small hump accommodating the positive skewness. For illustration purposes, we add in Fig. 2, both the non-Gaussian ME-pdf fit and the Gaussian pdf fit of the innovation histogram for the sample 4s. Asymptotically converging solutions to the ME-pdfs are obtained by solving the ME problem in finite interval sets D = [-L, L] of increasing length L up to stabilizing both the ME and the LMs. By using the Leibnitz derivation rule and fixing the expected values, it is easy to obtain the derivative of the ME H M with respect to L:

  a) sample size N ; b) Specified standard deviation σso, of the observation error, mean m(σ sb ) and standard deviation std(σ sb ) of the specified standard deviation of the background error, innovation bias b d , standard deviation σ d of the innovation (all in K); c) Skewness s d , kurtosis k d , and approximated negentropy J d of the innovation with the 99%-statistically significant values of s d ,k d and J d marked in bold; d) Heteroscedasticity measures h 1b , h 3b and h 4b , of σ sb ; e) Maximum variance fraction βg of the Gaussian error; f) Bounds of the tuning factors: f o1 , f b1 , f o2 and f b2 of the specified standard deviation of errors. See text for details concerning quantities in c), d), e) and f).

Figure 1 :

 1 Figure 1: Scatter plot of pairs: (specified background error standard deviation σ sb , North latitude) for channel 6. Boxes delimit samples corresponding to 'Ice-free' and 'Ice-covered' conditions.

Figure 2 : 30 - 3 Figure 3 :

 23033 Figure 2: Normalized histograms (in stairs), corresponding Gaussian pdf fit (grey solid line), and Maximum Entropy pdf fit (dark solid line) for the sample 4s of innovations.

Figure 4 : 32 - 4 Figure 5 :Figure 6 : 34 - 4 Figure 7 : 35 - 4 Figure 8 :Figure 9 :

 432456344735489 Figure 4: Contour map of the maximum variance fraction βg of a Gaussian error, as function of the innovation skewness s d and kurtosis k d . The contour map is symmetric for s d < 0

Figure 10 :Figure 11 :

 1011 Figure 10: Common to all figures a) to l): Solid thick line: Value of β b with the tuned mean variance of background errors. Three sets of curves with the scores MSC j (j = a, b, c) and non-Gaussianity source scenarios: (1) Gaussian background errors: black thin solid lines; (2) Gaussian observation errors: gray thin solid lines and (3) Maximum entropy sharing of non-Gaussianity: black dotted lines. The values of the parcial scores MSCa, MSC b , and MSCc for each scenario of non-Gaussianity and the same tuning factor f b obey in all cases to: MSCa < MSC b < MSCc. Figures a) to l) correspond respectively to samples 4s, 5s, 6s, 7s, 11s, 12s, 15s, 5i, 6i, 7i, 14i and 15i. The graphics are in logarithmic scale, values of the partial scores below 10 -4 are not plotted. See text for details.

Figure 12 :

 12 Figure 12: Lagrange Multipliers of the standardized ME-pdfs as function of its skewness s and kurtosis k: a) λ 1 ;b) λ 2 ; c) λ 3 ; d) λ 4 . e) Maximum entropy (nats) of the ME-pdfs in the SK domain, f) Standardized ME-pdfs for points G, SM, SP and SL.

  Table 1 lists the h-terms: h 1b , h 3b and h 4b for the 13 samples. The h 3b terms in Table 1 range in the interval [1.03, 1.29] from which one sees a more marked heteroscedasticity in samples 7s, 11s, 12s, 14s, and 15s. The deviations of h 1b and h 4b from unity follows that of h 3b among the available samples.


Jules Dupire

Manouk Abkarian

Annie Viallat

Chaotic Dynamics of Red Blood Cells in a Sinusoidal Flow

Keywords: numbers: 47.52.+j, 83.50-v, 83.80.Lz, 87.85.G time dependent flow, tanktreading, swinging, tumbling, capsule, elasticity, Lyapunov exponent

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Oscillators are ubiquitous in physics, chemistry (Belousov-Zhabotinski reaction) and in the living world at the scale of single organisms (pairs of beating eukaryotic flagella [START_REF] Goldstein | [END_REF]) or of populations (fireflies blinking). Ever since Huygens noticed in 1665 that two pendulum clocks synchronized, scientists discovered that coupled oscillators showed a range of complex behaviors from synchronous to chaotic oscillations. Recently, we demonstrated that the movement of a red blood cell (RBC) in a steady shear flow is well described by two ordinary equations analog to the problem of two dissipative coupled angular oscillators represented by the cell inclination angle relatively to the flow direction and the rotation angle of the membrane elements around the center of mass of the cell [2]. The two angles are synchronized, as required for 2D non-linear systems [2][3][4]. In the vasculature, however, the blood flow is time-dependent. In arteries it is pulsed and in the microcirculation, it presents irregular temporal variations, from intermittent cessation in capillary flows to chaotic behavior in arterioles [5][6][7]. These fluctuations, which are partly related to fluctuations of vascular calibre are still poorly understood. At the RBC scale, the time-dependent external flow plays the role of a third coupled oscillator, which, by adding a degree of freedom to the non-linear system could allow the cell to have a complex, even chaotic motion. Such motion would in turn influence the large scale flow of blood and contribute to its observed temporal fluctuations. In this work, we provide a theoretical description of the motion of RBCs in a sinusoidal moderate shear flow close to physiological shear stresses and flow variation frequencies that we illustrate with several experiments. We predict that RBCs can present a chaotic motion and we experimentally show the unstable movement of RBCs. We propose a full theoretical diagram of the RBC motion as a function of the amplitude and the frequency of the shear rate. We show that the existence of a cell shape memory associated with a minimum of the elastic energy of the membrane for a specific position of its elements is crucial for chaos. Finally, as our experiments allow an easy recording of RBC motion they may provide a suitable tool to characterize cell viscosity and elasticity.

Model -

The motion under a moderate shear flow of a RBC whose symmetry axis lies in the shear plane is characterized by two angles: θ, the cell inclination angle with respect to the flow direction and ω, the angle between the instantaneous position of a membrane element and its initial position at rest (Fig. 1). At low shear stress, RBCs tumble (T), i. e. θ rotates and ω oscillates. Above critical values of the external viscosity and the shear stress, the cell inclination stabilizes (θ oscillates around a mean value) while the membrane rotates around the fixed cell shape (ω rotates): the cell both tanktread (TT) and swing (S). The membrane elements, displaced during the tanktreading revolution are deformed and store an elastic energy which reaches a maximum for every π rotation of the membrane. This energy barrier modifies periodically the energy balance of the system and is at the origin the θ-oscillations (Fig. 1). At the T-STT transition, there is an intermittent regime with alternative series of T and STT [2,8,9], where the numbers of T and STT depend on the shear rate. A recent analytical simple model derived from Keller and Skalak's (KS) [10] was proposed [2,9], which treats a RBC as a viscoelastic ellipsoidal membrane enclosing a viscous fluid. It does not account for cell deformation and is therefore adapted to moderate shear stress. It deals with a basic mechanical description of the membrane and a simplified velocity field on the membrane, which lead to tractable analytical equations that retrieved the observed regimes of cell motion and provided semi-quantitative fits of experimental data [2]. Subsequent full numerical models [11][12][13] derived from this approach basically retrieved the same features for the cell motion. The RBC shape is described by a fixed oblate ellipsoid filled with a Newtonian solution of viscosity η i and delimited by a viscoelastic thin membrane, which obeys a simple Kelvin-Voigt constitutive law: σ = 2η m D + 2µ m E, where η m and µ m are the membrane viscosity and shear modulus and D and E are the strain rate and strain tensors, respectively. The sinusoidal shear rate of the flow is γ = γa sin(νt). The balance of torques exerted by the external fluid on the cell and the conservation of energy [10] -the rate of viscous dissipation in the cell plus the elastic power stored by the cytoskeleton is equal to the power provided by the external fluid on the cell -yield the two equations for RBC motion given in ref [2]. These equations still hold for time-dependent shear rate. Adding a third equation with ν = νt yields a 3D autonomous flow F ( X) such as d X/dt = F ( X), where X is the vector (θ, ω, ν): For small shear rate ( γa < γc ), where γc is the critical shear rate of the steady T-STT transition, RBCs tumbles periodically. The phase trajectories approach a stable limit-cycle as shown in Fig. 1a on the projection on the plane (θ, ω). The approach to the limit-cycle takes about 1 sec. The angle θ changes monotonously with time: clockwise during the forward movement ( γ > 0) and counterclockwise back to its initial value during the backward movement ( γ < 0). The angle ω oscillates, indicating that the membrane slightly moves back and forth about its position at rest.

∂ t θ = γa sin(ν) - 1 2 + a 2 1 -a 2 2 a 2 1 + a 2 2 cos(2θ) - 2a 1 a 2 a 2 1 + a 2 2 ∂ t ω ∂ t ω = f 3 γa f 2 -ηi ηo (1 + ηm ηi Ω V )f 1 (1) × f 1 2f 3 µ m η o γa Ω V sin(2ω) -sin(ν)cos(2θ) ∂ t ν = ν f 1 ,
For larger values of the shear rate amplitude ( γa > γc ), the cells swing and tanktread when the instantaneous shear rate is larger than γc and tumble when it decreases below γc just before the flow changes direction, as shown on an experimental sequence in Fig. 2. The model predicts a stable periodic regime described by either one or two non-intersecting limit-cycles (see Fig. 1b), whose trajectories are determined by the initial angles. During the S-TT regime, θ is little sensitive to the initial conditions. However, the storage of elastic energy slightly differs on each trajectory. Therefore the TT-T transition does not occur at the same moment and the values of θ and ω when the flow changes direction differ from one trajectory to the other. The consequence for the cell reversal is major: in one case, at the end of the forward movement the cell has a counterclockwise reversal and passes vertically (see Fig. 2), perpendicularly to the flow direction, whereas in the other case, the cell rotates in the clockwise direction and passes horizontally, in the alignment of the flow direction.

For specific values of γa and ν, a non-periodic motion, highly sensitive to initial conditions and characterized by an unstable limit-cycle occurs. Projections of portions of trajectories are shown in Fig. 1c. Trajectories are no longer closed and two initially very close initial conditions lead to strongly divergent trajectories. For illustrations of the attractor and first return map see [14]. In this case, the three oscillators θ, ω and γ cannot synchronize and the movement is chaotic. Chaos was further investigated by computing the Lyapunov exponents of the system, which characterize the rate of separation of initially close trajectories. Let X(t) be a solution of the flow F and δ X(0) an initial infinitesimal separation. The separation at time t, δ X(t), obeys the equation of evolution obtained by linearisation of the flow F , and its integration leads to the exponential evolution of the perturbation δ X(t):

d δ X dt = ∂F ∂ X X(t) δ X; δ X(t) = δ X(0) e R t 0 ∂F ∂ X | X(t) dt (2)
For each given solution, we computed the Jacobian matrix ∂F/∂ X, we numerically determined the matrix

t 0 ∂F ∂ X | X(t)
dt, and we calculated its three eigenvalues at several times. The slopes of the asymptotic linear variations of these values with time are the three Lyapunov exponents of the system. The one corresponding to the parameter ν is equal to 0 and the second one is always negative. The largest exponent λ may be positive (see [14]), which is a signature of chaotic motion. Variations of λ and of the error on λ (determined from the correlation coefficient of the least square fit, 1 -r 2 ) with the shear rate frequency are shown in the insert of Fig. 3. The error is small for the largest positive values of λ, thus confirming the existence of chaotic zones of the RBC motion. We analyzed the Floquet matrix for the Poincare section computed when the instantaneous shear rate equals 0 at the bifurcation event. The Floquet eigenvalue crosses the unit circle along the real axis at + 1, wich characterizes a saddle-node bifurcation and discloses that the switch to chaotic behavior occurs via an intermittency (type I) mechanism, characterized by bursts of chaos. We fully characterized the RBC dynamical behavior by plotting λ as a function of ν and the RMS value of the shear rate γrms = γa / √ 2 (Fig. 3). Values for ν span physiological frequencies, from vasomotion fluctuations (tenths of Hertz) to the heart frequency (one Hertz) and values of γ cover the range of stresses η o γ observed in the vaculature, from capillaries to arteries [5][6][7]. When γrms decreases below 3.1, λ deeply falls around -0.5, revealing that the system is stable. The value 3.1 is equal to γc , the critical shear rate of the T-TT transition in constant shear flow for the RBC we computed. Tumbling stabilizes the motion. Indeed, during tumbling, the cell membrane is 'quasi-solidified', the angle ω remains very small and it can be eliminated in the first equation of eq. 1 by using the second equation where the term in sin(2ω) is neglected. The system then loses one dimension and chaos is no longer possible. Above γc the motion is much less stable. Unstable chaotic zones appear. They form widening bands located on a bundle of lines. The solid line represents the ( γrms , ν) couples for which a RBC performs a full swinging cycle during the tanktreading time of the forward motion ( γ > 0) and one swinging cycle during the backward motion ( γ < 0). The other dotted lines characterize the flow conditions for which a RBC swings respectively 1/2 (steepest line), 2, 3 and 4 times during half a flow period. The strong correlation between these lines and the chaotic bands shows that chaos appears at the resonance of the flow frequency with the internal specific cell frequency. As expected, chaos is intimately related to swinging, i. e. to the elastic storage of the non spherical membrane, since the term with the elastic modulus µ m of Eq.1 induces the coupling of the three differential equations. Viscous capsules or spherical viscoelastic capsules must present a stable dynamics in a time-dependent shear flow.

Observations -Direct measurements of θ are provided from side-view microscopic imaging in a vertical plane parallel to the shear plane [2,15]. Sinusoidal flow was created from the variation of hydrostatic pressure by connecting the flow chamber to a reservoir (10ml) bound to the outer edge of a vertical rotating wheel (radius 9.45 cm). The frequency of rotation of the wheel ν (in rd/s) was varied in the range [0.16 -1.15rd/s] and γa = 6.72 s -1 . RBCs were diluted in a solution of dextran (concentration 9% w/w, η o = 34 mPa.s, Mw= 2.10 6 g/Mol) with PBS at 290 mOsm and pH=7.4. RBCs were almost non-buoyant in the solution and could flow during several hours without sinking. The membrane shear modulus as well as the internal and membrane viscosities of RBCs depend on the age of the cell and are therefore widely distributed when measured in a blood sample. The maximal applied shear stress η o γa was specifically chosen to be larger than the critical shear stress η o γc for the majority of RBCs, so that most cells swing and tanktread but the stiffest ones tumble. We recorded cell movements 10 minutes after application of the flow in order to let the motion stabilize and reach the limit-cycle. We observed the three motions predicted by the model. Stable tumbling is characterized by the continuous variation of θ with time as illustrated in Fig. 4a over three flow 4 and in [14] for a single cell observed at 6 different flow frequencies. In the case of unstable behavior, an adjustement is stricto sensu not possible since the motion is highly sensitive to initial conditions but the theoretical curve well describes the main characteristics of the motion (see inserts of Fig. 4) and predicts the instability of motion.

Here, we predict that dilute RBCs present a chaotic motion over large domains of external sinusoidal flows, under shear stresses and flow variation frequencies relevant for arterial and microcirculation physiology. In-vivo, in the whole blood cells form a concentrated suspension and interact together. These interactions couple more oscillators together and should result in enhanced chaotic motions which could affect flow properties in the sheared blood layer close to vessel walls. Finally, our approach also applies to other physical systems used in flow as polymer capsules or protein-coated drops.

f 2 and f 3

 3 are geometrical constants described in ref[10], a 1 and a 2 are the lengths of the long and small axis of the cell cross-section, V is the RBC volume, Ω is the membrane volume and η o is the external suspending fluid viscosity. The phase space is conveniently described by the ensemble (θ, ω, sin(ν)). Equations are solved numerically by using the set of parameters[2]: a 1 = a 3 = 4µm, a 2 = 1.5µm, Ω V = 7.48.10 -2 , η 0 = 34.10 -3 Pa.s, η i = 10 -2 Pa.s, η m = 0.7 Pa.s, µ m = 1.6 Pa.

FIG. 1

 1 FIG.1: a-c: numerical simulation of θ( • ) vs. ω( • ). a: limit cycle γa = 2s -1 , ν = 0.2 rd.s -1 . b: two limit cycles γa = 9s -1 , ν = 0.35 rd.s -1 ; Initial conditions: ω0 = 0, θ0 = 0 (-) or θ0 = π/2 (--). c: non-closed trajectory for γa = 9s -1 , ν = 0.8 rd.s -1 , initial conditions: ω0 = 0, θ0 = 10 -9 (-) or θ0 = -10 -9 (--). d: sketch of a RBC in the flow. e: numerical chaotic trajectories θ( • ) vs. time, initial conditions of c.

FIG. 3 :

 3 FIG. 3: Diagram of the dynamical RBC motion; a: 3D map of λ vs. ν and γrms; b: 2D projection: values of λ (in color) in the (ν, γrms) plane. Dark regions (positive λ) are chaotic zones. Solid line: flow conditions for which RBCs perform a full swinging cycle during the forward (or backward) motion. Dotted lines from top to bottom: RBCs perform respectively 1/2, 2, 3 and 4 swinging cycles during the forward motion; insert: variation of λ and 1 -r 2 with ν in a chaotic zone.

FIG. 4 :

 4 FIG. 4: Experimental variation and fits of θ vs. time for 3 RBCs. a: ν = 0.08rd.s -1 , stable T, fit η ef f = 1.25 10 -2 P a.s, µ = 15.2P a; b: ν = 0.05rd.s -1 , stable S-TT and T, fit η ef f = 4.95 10 -2 P a.s, µ = 0.591P a; c: ν = 0.11rd.s -1 . Unstable chaotic motion, fit η ef f = 3.97 10 -2 P a.s, µ = 0.49P a. Inserts show a zoom of the boxed regions

Struct.,23, 787 (1994).

We thank Dr. R. Bruinsma for his suggestions, P. Soulé for discussions and C. Meuzard for experiments. Our labs belong the CNRS French Consortium 'CellTiss'.