Marianna Meo
email: marianna.meo@ihu-liryc.fr

Olivier Meste
email: meste@i3s.unice.fr

Sergio Signore
email: sergiosignore@gmail.com

Marcello Rota
email: marcello_rota@nymc.edu

Novel Methods for High-resolution Assessment of Cardiac Action Potential Repolarization

Keywords: Diabetes, murine model, singular value decomposition, action potential, cardiac repolarization

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

Introduction

The action potential (AP) of cardiomyocytes underlies the electrical and contractile behavior of the entire heart, with important implications on the stability of cardiac rhythm and mechanical function. The analysis of the AP of myocytes provides information on the electrical behavior of cardiac cells obtained from the normal or diseased heart, and modalities of adaptation of the electrical activity to experimental interventions. Thus, the examination of the AP of myocytes under various circumstances allows the identification of alterations contributing to the origin of electrical and mechanical disturbances, on the one hand, and provides information on putative ionic currents underlying changes of electrical properties of cardiac cells, on the other.

Phase-plane analysis of cardiac AP has been previously proposed [START_REF] Sperelakis | Phase-plane analysis of cardiac action potentials[END_REF] to graphically compute depolarization rate constants and equilibrium potentials, and to obtain information on ionic currents responsible for the profile of the AP. However, the analysis of phase-plane plots is merely visual, and the quantification of properties of AP duration and morphology is geometrically and computationally complex, limiting the applicability of this analysis when evaluating APs obtained from different myocytes or from the same cell under various conditions.

Currently, APs are analyzed using semiautomatic programs allowing the computation of the duration of the AP (APD) at discrete, user-defined repolarization levels, while neglecting information present in other ranges of the electrical recovery. This aspect reduces the level of accuracy of the assessment and questions the suitability of currently available approaches to capture temporal adaptations and beat-to-beat variations of the repolarization phase of the AP.

To address these limitations, we have implemented a novel method for the assessment of repolarization properties of the AP and alterations occurring during a train of stimulations. Initially, the inverse function of the electrical signal of the AP was computed in order to express duration of the AP as a function of the repolarization level. Subsequently, alterations occurring in the repolarization phase were analyzed using two distinct approaches. With the first approach, beat-to-beat alterations of the AP repolarization were either computed directly on the raw data or determined through linear or higher order polynomial regressions of APD values. With the second approach, a modified evolutive singular value decomposition (ESVD) method was applied to APs, and information related to their temporal dynamics (i.e., changes occurring in consecutive APs) was obtained for the entire repolarization range.

To test the functionality of our methodologies, AP repolarization properties were assessed in myocytes and perfused hearts from healthy and diabetic mice.

This selection was based on previous finding documenting that hyperglycemia is associated with prolongation of the AP and enhanced dispersion of cardiac electrical repolarization [START_REF] Meo | Reduction in Kv current enhances the temporal dispersion of the action potential in diabetic myocytes: insights from a novel repolarization algorithm[END_REF][START_REF] Sorrentino | Hyperglycemia induces defective Ca2+ homeostasis in cardiomyocytes[END_REF], possibly interfering with temporal adaptations of the AP profile. We report that the modality of adaptation of the AP repolarization to a train of stimulations differs in control and diabetic myocytes. Thus, our novel approaches represent valuable platforms for the detailed and comprehensive assessment of AP properties of cardiomyocytes.

Materials and Methods

Experimental animals, patch-clamp and ex-vivo studies

Female C57Bl/6 mice at 3 to 11 months of age were studied in accordance with the Guide for Care and Use of Laboratory Animals; experiments were approved by the local animal care committee. Diabetes was induced by streptozotocin (STZ, Sigma-Aldrich, St. Louis, MO) injection, as previously reported [START_REF] Meo | Reduction in Kv current enhances the temporal dispersion of the action potential in diabetic myocytes: insights from a novel repolarization algorithm[END_REF][START_REF] Sorrentino | Hyperglycemia induces defective Ca2+ homeostasis in cardiomyocytes[END_REF]. Animals with blood glucose level higher than 300 mg/dL were included in the study. Data were obtained from cardiomyocytes collected from hearts of naïve mice (control, CTRL, n = 41 cells from 14 animals) and STZtreated hyperglycemic mice (STZ, n = 76 cells from 10 animals) [START_REF] Meste | A new tool for the action potential repolarization dynamic analysis: application to the discrimination of diabetic and control cells[END_REF]. Details related to the enzymatic digestion of left ventricular myocytes from murine hearts and current-clamp experimental setup are provided elsewhere [START_REF] Meo | Reduction in Kv current enhances the temporal dispersion of the action potential in diabetic myocytes: insights from a novel repolarization algorithm[END_REF][START_REF] Signore | Late Na(+) current and protracted electrical recovery are critical determinants of the aging myopathy[END_REF][START_REF] Borghetti | Notch signaling modulates the electrical behavior of cardiomyocytes[END_REF].

Monophasic action potentials (MAPs) were recorded at a 4 kHz sampling rate using a micro MAP-Tip electrode (Harvard Apparatus) positioned on the left ventricular free wall of Langendorff perfused hearts from 7 CTRL and 5 STZ mice. Hearts were stimulated via a mini-coaxial electrode at 125-ms basic cycle length. Details related to the experimental setup are provided elsewhere [START_REF] Meo | Reduction in Kv current enhances the temporal dispersion of the action potential in diabetic myocytes: insights from a novel repolarization algorithm[END_REF][START_REF] Signore | Late Na(+) current and protracted electrical recovery are critical determinants of the aging myopathy[END_REF].

Formatting of electrical recordings

For our analysis, a series of 70 consecutive APs obtained from CTRL and STZ myocytes stimulated at 2 Hz pacing rate was employed. Electrical signals were initially processed using Clampfit 10 software (Molecular Devices) for calculation of the duration of the AP at 10% (APD10), 30% (APD30), 50% (APD50), 70% (APD70), and 90% (APD90) of the repolarization phase. In brief, APDR was computed as the time between the positive peak of the AP and the level when the AP dropped (repolarized) by R% of its height. Electrical signals were converted into text files for processing in MATLAB (MathWorks).

APs were segmented and aligned with respect to the depolarization time and arranged into a I × N matrix X. I = 70 signifies the number of APs examined in a recording and N = 5000 reflects the number of samples in each stimulation cycle, corresponding to 500 ms:

X = [x 1 , x 2 , . . . , x I ] T , (1)
with each vector x i , i = 1, . . . , I representing a single AP. Representative APs obtained from a CTRL and STZ myocytes are reported in Fig. 1 (upper panel).

Data were acquired at a sampling rate of 10 kHz and, subsequently, electrical recordings were downsampled by a factor of 10, producing N * = N/10 samples. This step was performed to reduce the computational load without loss of relevant information.

Analysis of temporal alterations of the repolarization phase of the AP

To quantify alterations in the repolarization phase of the AP occurring during a train of stimulations and to compare the behavior observed in CTRL and STZ myocytes, the following steps were applied (Fig. 2).

Computation of the inverse function of AP signals

To present duration of the AP as function of repolarization, AP signal matrix was processed as reported previously [START_REF] Meo | Reduction in Kv current enhances the temporal dispersion of the action potential in diabetic myocytes: insights from a novel repolarization algorithm[END_REF][START_REF] Meste | A new tool for the action potential repolarization dynamic analysis: application to the discrimination of diabetic and control cells[END_REF][START_REF] Meste | Diabetes affects the temporal dynamics of the repolarization properties of cardiomyocytes[END_REF] (see Fig. 1, lower panel).

AP amplitude was normalized in the repolarization range [0% -95%] between 1 and 0. The terminal portion of the AP, [95% -100%], was not processed to limit spurious artifacts due to boundary conditions. After this transformation, the late repolarization appeared on the left portion of the horizontal axis, whereas the early phase of repolarization appeared on the right side (see Fig. 1, lower panel). After this signal conditioning, the notation X was maintained for the observations' matrix. Computing the inverse of x i required the function to be strictly monotonic, a condition which is met by the rodent AP.

Each AP observation x i , i = 1, . . . , I was modeled as a function of its duration n, n = 1, . . . , N * , as:

x i (n) = f i (n;β β β i ) + e i (n) (2)
where f i (n;β β β i ) is a parametric function with parameters' vector β β β i , and the term e i (n) accounted for noise. As previously done [START_REF] Meste | Diabetes affects the temporal dynamics of the repolarization properties of cardiomyocytes[END_REF][START_REF] Cabasson | Quantifying the PR interval pattern during dynamic exercise and recovery[END_REF], f i (n;β β β i ) was expressed as a linear combination of a collection of functions w l (n):

f i (n;β β β i ) = L l=1 β i,l w l (n) (3)
where the number of functions L was assumed as the maximum integer not greater than (N * + 1)/2 (see Supplementary Material in [START_REF] Meo | Reduction in Kv current enhances the temporal dispersion of the action potential in diabetic myocytes: insights from a novel repolarization algorithm[END_REF] for graphical explanation). As previously reported [START_REF] Meste | Diabetes affects the temporal dynamics of the repolarization properties of cardiomyocytes[END_REF][START_REF] Cabasson | Quantifying the PR interval pattern during dynamic exercise and recovery[END_REF], functions w l (n) were chosen as piecewise linear with a triangle shape. It was demonstrated that imposing f i (n;β β β i ) to be strictly monotonically decreasing implied that the same constraint could be extended to β β β i as well:

∀l ∈ [1 : L -1], β i,l > β i,l+1 > 0 (4)
By using this condition, β β β i could be determined by least squares estimation under inequality constraints, with conversion to a least distance programming (LDP) problem [START_REF] Lawson | Solving least squares problems[END_REF]:

β β β i = arg min β β β i N n=1 (x i (n) -f i (n,β β β i )) 2 (5)
Subsequently, AP observations x i were replaced by the new estimated function f i , which was strictly decreasing, and it appeared smoothed, because the noise e i had been filtered out. By this process, the inverse function n describing the profile of the AP was uniquely interpolated as:

n = f -1 (z,β β β i ), ∀z ∈ [0, 1]. (6)
In this study, M = 500 repolarization level values z were uniformly selected in the interval [0, 1], providing sufficient accuracy for subsequent processing. These inverted APs were then analyzed by the two newly developed approaches.

Mathematical models of AP repolarization profile

Following computation of the inverse function for electrical signals, raw AP duration data were employed to evaluated temporal alterations occurring during the stimulation train.

Subsequently, we investigated whether it was possible to determine a more general mathematical function which could accurately model the duration of AP (APD) computed in Eq. ( 6) over the entire repolarization range and across the course of the stimulation. We first explored the suitability of a linear law to quantify the relation APD vs stimulation index, for each repolarization level.

Then, we investigated whether the inclusion of additional, non-linear terms could improve the descriptive power of the model and capture more complex phenomena in the cardiac electrical recovery. For this purpose, we tested several polynomial functions with increasing regression orders.

The use of linear and higher order polynomial functions to describe the duration of AP over the repolarization range was selected because of the possibility to quantify APD changes in response to the stimulation in terms of their derivative, which can be easily computed and expressed in terms of the polynomials' coefficients. More precisely, these coefficients can provide information on the direction of APD temporal variations, i.e., on the overall shortening/lengthening of repolarization along the stimulation, allowing a more straightforward and clear interpretation of the data.

The strategies presented so far aim to determine APD changes over the progression of the stimulation by considering each repolarization level separately. As an alternative approach, we determined a common behavior of the repolarization phase of the AP over the stimulation train and then quantified the deviation from this common pattern, for each repolarization level r. To this end, we investigated the applicability of ESVD, which allows the extraction of the most meaningful signal components, providing a compressed representation of input data and easing their interpretation. With regard to the analysis of AP electrical recovery, by ESVD we extracted an AP template summarizing the main trend of repolarization (i.e., overall APD increasing/decreasing), and then used the obtained template as a reference for a more detailed analysis of changes of the AP pattern during the stimulation, at specific repolarization levels. Compared to standard SVD, ESVD can provide insights into subtle temporal variations of AP observations, which can be easily quantified in terms of function derivatives.

In the next subsections, more mathematical details about these AP models are provided.

Assessment of temporal trends of AP repolarization by direct data processing

To evaluate alterations of the AP occurring during the stimulation train, durations of AP (APDs), computed with Eq. ( 6) at the r-th fixed level of repolarization r, were either directly processed or interpolated via linear or higher order polynomial regression. At each repolarization level r = 1, . . . , M , a set of I = 70 APD values (one per beat) was measured, and each sequence was denoted as a vector y r , with r = 1, . . . , M .

Values of the AP duration over the course of the stimulation for a given repolarization level were fitted with linear regression or polynomial functions up to the fifth order. Trends in the duration of the AP were quantified by the mean of the first derivative of raw data or relations obtained by linear or higher order polynomial function. The mean of derivative values was then normalized by the mean of the APD for that repolarization level and labelled as APT DIR .

APT DIR was employed as an index of directionality and rate of APD change.

Assessment of temporal trends of AP repolarization by ESVD

As previously reported [START_REF] Meste | A new tool for the action potential repolarization dynamic analysis: application to the discrimination of diabetic and control cells[END_REF], we introduced the first order model for r = 1, . . . , M repolarization levels:

y r (i) = p r (i).v(i) + e r (i), (7)
where i = 1, . . . , I and e r (i) accounts for the acquisition and noise modeling.

The problem presented here is an order-one approximation of the more complex decomposition:

y r = p 1,r • v 1 + p 2,r • v 2 + p 3,r • v 3 + . . . ( 8
)
which, in turns, is an extension of the standard SVD:

y r = a 1,r v 1 + a 2,r v 2 + a 1,r v 3 + . . . ( 9
)
This approach is based on the functional data analysis theory [START_REF] Madrigal | Uncovering correlated variability in epigenomic datasets using the Karhunen-Loeve transform[END_REF], assuming that measurements originate from continuous phenomena. The function p r (i) is related to the r-th repolarization percentage, whereas the term v(i) contains information common to all levels r = 1, . . . , M . Functions p r (i) were decomposed over a set of K basis functions b k (i) such that:

p r (i) = K-1 k=0 b k (i)θ r,k , (10)
In the sequel, a set of K=3 polynomial functions b k (i) was introduced in order to obtain smooth variations of the AP profile when varying i over the repolarization interval. In vector format with samples arranged column-wise, Eq. ( 7) can be expressed as:

y r = p r • v + e r = B • (v1I T )θ θ θ r + e r = M v θ θ θ r + e r ( 11
)
where the coefficients b k (i) are stored column-wise in B, θ θ θ r contains the θ r,k values, and As previously performed [START_REF] Meste | A new tool for the action potential repolarization dynamic analysis: application to the discrimination of diabetic and control cells[END_REF], polynomials p r are modeled as 3rd-order functions and considered multiplicative factors of vector v, thus dimensionless. To assess AP profile variations for each repolarization level r in the range [0% -95%], the mean of the first derivative of p r over the course of the stimulation train was determined, together with the sign of the mean of the first derivative of v.

M v = B • (v1I T ).
The first term refers to AP dynamics specific to a fixed level of repolarization r, whereas the second term describes an average behavior which is common to all phases of repolarization. The product of the two, denoted APT ESVD , exhibits the following property: if APT ESVD >0 then p r reinforces the trend of v, i.e., the extent of the AP profile to shorten/prolong along the stimulation.

By contrast, if APT ESVD <0 the term p r weakens the trend of v.

Assessment of cardiac repolarization at the tissue level through automatic MAP processing

To determine whether the analysis of the electrical recovery at the single cell level through the two new algorithms could be extended to the intact myocardium, MAPs were recorded in Langendorff perfused hearts obtained from n=7 CTRL and n=5 STZ mice. A set of 40 consecutive beats from these electrical signals was processed. The average MAP profiles were determined for both groups. Furthermore, APT was assessed by both methodologies, i.e., the direct approach (with and without polynomial interpolation) and the ESVD approach.

Statistical analysis

Data were initially tested for normality with Agostino-Pearson test for assignment to parametric or non-parametric analysis. Accordingly, significant deviations from zero were verified through t-test or Wilcoxon signed rank test.

Similarly, temporal changes in the AP indices throughout the stimulation train were assessed via paired t-test or Wilcoxon matched-pairs signed rank. For each method, i.e., standard APD analysis, direct processing (with or without polynomial approximations), and ESVD processing, intergroup differences were evaluated via rank-sum test. A similar approach was applied to repolarization indices computed from MAP recordings. Statistical tests were considered significant for p<0.05.

Results

Temporal alterations of the AP repolarization during a train of stimulations

To test the ability of our algorithms to define temporal alterations occurring in the repolarization phase of the AP and identify putative differences in the behavior of cardiomyocytes from normal and diseased hearts, 70 sequential APs elicited during a train of stimulations applied to CTRL and diabetic cells were studied. Specifically, the sequence of APs was obtained in 41 CTRL and 76 STZ myocytes. Initially, we evaluated whether changes in AP properties occurred during the sequence of stimuli. Therefore, duration of the AP at 10% (APD10), 30% (APD30), 50% (APD50), 70% (APD70), and 90% (APD90) of the repolarization phase was computed for the first and last AP of the stimulation sequence. Analysis was computed for each cell of the CTRL and STZ groups using standard methodology (Clampfit software) and computation of the inverse AP function, as previously done [START_REF] Meo | Reduction in Kv current enhances the temporal dispersion of the action potential in diabetic myocytes: insights from a novel repolarization algorithm[END_REF] (Fig. 3). The two methodologies yielded to similar results, as documented by linear regression and Pearson correlation analysis for the two datasets (Fig. 3A and Fig. 3B). With data obtained with the two approaches, durations of the AP at different repolarization levels for the first and last beats were then compared. Quantitatively, for CTRL myocytes, APD10 was significantly prolonged at the end of the stimulation train, with respect to the initial beat (Fig. 3C). Similarly, for diabetic cells, APD10 and APD30 were prolonged at the end of the stimulation train, but shortening of the AP occurred at 70% and 90% of the repolarization phase (Fig. 3D). Thus, based on observed temporal changes in AP properties for CTRL and STZ myocytes, electrical signals were processed using our newly developed analytical approaches to determine whether these methodologies offered additional information.

Evaluation of AP repolarization properties by direct processing of electrical signals

To establish the ability of our algorithms to define temporal alterations occurring in the repolarization phase of the AP and identify putative differences in the behavior of CTRL and diabetic cardiomyocytes, APs elicited by the train of stimulations in the two groups of cells were processed. Initially, we compared outcomes of the processing stage involving linear regression or higher order polynomial interpolation of APD over the course of the stimulation train, with respect to the raw input data. A representative example of different approximations for the duration of the AP at 55% and 95% of the repolarization for a STZ myocyte over the course of a train of 70 stimuli is reported in Fig. 4. For these two repolarization levels, the AP tended to become shorter with progression of the stimulation train (stimulation index). Linear regression and or higher order polynomial interpolation of APD input raw data successfully attenuated beat- to-beat variability caused by data structure and acquisition noise. Importantly, temporal trends of AP repolarization at 55% and 95% were preserved and properly characterized by the two regression approaches. With respect to the linear interpolation, increasing the regression order up to 5 did not bring any substantial change to the characterization of APD behavior. The similarity between the raw profiles and the regression curves for APD at the two repolarization levels corroborates the effectiveness of the proposed approach to describe alterations occurring in the AP over the course of a stimulation train. Subsequently, we evaluated whether the direct processing stage allowed a quan- tification of the response of CTRL and STZ cells to the train of stimulation (temporal trends). For this purpose, the first derivative for regression curves of APD was obtained for the entire repolarization phase; their mean was calculated and normalized by the mean duration of the APD (APT DIR ). This processing stage allowed us to define direction and magnitude of AP changes over time. Overall, for both CTRL and STZ myocytes, APT DIR differed from zero, indicating that the repolarization phase was altered (i.e., became longer or shorter) during the stimulation train (Fig. 5, upper panel). Then, signed rank test median statistic of APT DIR was employed to establish whether detected temporal trends for the repolarization profile of CTRL and STZ myocytes were statistically significant (Fig. 5, lower panel). For both group of cells, APT DIR was positive in the early repolarization phase, consistent with a progressive prolongation of this portion of the AP. In contrast, APT DIR was negative in the late repolarization, reflecting a progressive shortening of this phase of the AP for CTRL and STZ myocytes. Importantly, the transition of APT DIR from positive to negative values occurred at ∼ 40% and ∼ 60% of the repolarization phase for the AP of CTRL and STZ cells, respectively, indicating that adaptations of the AP profile to a stimulation train tended to be different in the two groups of cells. However, statistical analysis by unpaired Wilcoxon rank-sum test of APT DIR plots computed from raw data, linear regression, or 5th-order polynomial interpolation revealed that differences in the temporal adaptation of the AP of CTRL and STZ myocytes did not reach statistical significance (Fig. 6). Overall, our methodology involving direct processing of electrical signals allow a comprehensive quantification of alterations occurring in the profile of the AP during a train of stimulation.

Evaluation of AP repolarization properties by ESVD of electrical signals

The use of matrices for the computation of the AP profile over the course of stimulation allowed us to implement matrix decomposition approaches to quantify changes of the AP. Thus, a signal processing involving evolutive singular value decomposition (ESVD) method was employed to establish the ability of this approach to capture temporal adaptations occurring in CTRL and STZ myocytes. In Fig. 7 By rank-sum test of APT ESVD performed on the entire repolarization, it was found that, by ESDV processing, the temporal adaptation of the AP of CTRL and STZ myocytes, was significantly different only in the interval ∼50 -65% of the repolarization phase (Fig. 9). Thus, ESDV analysis provides complementary information with respect to direct processing of electrical signals.

Repolarization algorithms and electrical signals collected in the intact myocardium

To establish whether processing of electrical signals by direct approach and ESVD methodologies were applicable to extracellular recordings obtained in the intact myocardium, a sequence of monophasic action potentials (MAPs) recorded in Langendorff perfused hearts from CTRL and STZ mice were analyzed. The mean profile of the MAP repolarization was prolonged in STZ hearts (n=5), with respect to CTRL (n=7) (Fig. 10, upper panel). Consistent with previous observation [START_REF] Meo | Reduction in Kv current enhances the temporal dispersion of the action potential in diabetic myocytes: insights from a novel repolarization algorithm[END_REF], the difference was statistically significant in the range of 15%-55% of the repolarization phase (Fig. 10, lower panel). Temporal alterations for the repolarization phase of MAPs were then assessed in 40 consecutive beats for hearts stimulated at 8Hz, by direct and ESVD approaches. For each heart, we elected to restrict our analysis to 40 MAPs to avoid inclusion of electrical signals contaminated by motion artifacts and/or probe displacement occurring during the acquisition. The direct approach revealed a tendency for a prolongation of the early phase of the MAP in CTRL and STZ hearts, whereas the late repolarization remained stable, or slightly shortened (Fig. 11,A). The prolongation of MAPs for CTRL hearts reached statistical significance in a small portion of the early repolarization phase (Fig. 11,C). Overall, no statistical differences were observed between CTRL and STZ hearts for alterations occurring in MAP profile during the course of stimulation (data not shown). In a comparable manner, ESVD approach revealed alterations in the early and late repolarization phases of MAPs for CTRL and STZ heart (Fig. 11,B), but these changes were not significant (Fig. 11,D). Importantly, statistical comparison indicated that temporal behavior of MAPs for CTRL and STZ hearts differed, collectively, for ∼ 10% of the repolarization phase. These findings reflect, at least in part, observations at the cellular levels.

Thus, new developed methodologies are appropriate to process electrical signals collected in isolated cardiomyocytes and intact tissue. we graphically displayed repolarization levels on the horizontal axis, and determined the corresponding duration on the vertical axis. The entire range of AP repolarization was analyzed with high resolution, and was not restricted to a limited set of user-defined, discrete levels, as commonly done. Moreover, the inversion preprocessing stage enabled the identification of patterns related to the duration of the repolarization phase and beat-to-beat variations. This is in contrast with time-domain approaches presenting transmembrane potential as a function of time, an aspect that renders these methodologies ideal for the analysis of amplitude and rates of depolarization/repolarization [START_REF] Hodgkin | The effect of temperature on the electrical activity of the giant axon of the squid[END_REF] rather than direct APD quantification.

Previously, we demonstrated that mean duration and variability of AP could be effectively determined through the preprocessing stage involving the inversion of the relation between transmembrane potential and time [START_REF] Meo | Reduction in Kv current enhances the temporal dispersion of the action potential in diabetic myocytes: insights from a novel repolarization algorithm[END_REF]. Specifically, our algorithm yielded results comparable to those obtained via standard approaches when evaluating APD at discrete levels of the repolarization phase. Results obtained here further corroborate the ability of our method to properly report AP repolarization characteristics and their modifications over time. The level of detail and resolution of our approach cannot be achieved with standard methodologies. With respect to the assessment of temporal alteration of repolarization, our direct and ESVD approaches accurately described alteration in AP profile and characterized the impact of a pathological condition on AP repolarization.

Specifically, the direct approach identified that in control and diabetic cells APD tended to increase in the early phase of repolarization, whereas it pro-gressively decreased in the late phase. However, the transition from progressive lengthening to progressive shortening of the duration of the AP was identified at different levels of the electrical recovery for CTRL and STZ myocytes. These results are consistent, at least in part, with physiological findings presented by Nobe and colleagues [START_REF] Nobe | Chronic diabetes mellitus prolongs action potential duration of rat ventricular muscles: circumstantial evidence for impaired Ca2+ channel[END_REF]. They reported that transition in APD dynamics occurred between 50% and 75% of the repolarization in rat cardiomyocytes, in the presence of the inhibition of calcium release from the sarcoplasmic reticulum.

However, due to the limited resolution of the analysis, a precise evaluation of the transition level was not possible. In contrast, our high-resolution method captured changes in APD variation trend at 60% and 40% of repolarization for cells from STZ and healthy animals, respectively. Additionally, ESVD analysis revealed that, at the cellular level and in the intact myocardium, temporal behavior of the intermediate/late repolarization phase was different in CTRL and STZ groups. This feature would have not been identified when considering discrete, user-defined repolarization levels. Overall, our proposed methodologies offer higher detail in the assessment of AP repolarization properties with respect to commonly used, standard approaches.

Study limitations and future implementations

With this investigation we have documented the ability of our methodologies to detect temporal alterations occurring in the repolarization phase of the AP of isolated cardiomyocytes. Also, we have provided initial results on the suitability of our repolarization algorithm to characterize the profile of MAP and to monitor putative temporal alterations of the electrical activity of the intact myocardium. Thus, our results, which are merely methodological in nature, highlight the possibility to employ this computational platform for the analysis of electrical recordings collected in single cells, multicellular preparations, and myocardium.

Isolated cells lack the influence of cell-to-cell interactions and electrical coupling present in the myocardium. Thus, whether reported temporal alterations of AP profile observed in isolated cells are present and occur with a comparable magnitude in vivo and in the intact myocardium remains to be properly addressed. We found minor temporal alterations for MAPs recorded in the entire heart. However, differences in pacing paradigms for our studies at the cellular and organ levels preclude a direct comparison of these two datasets. The effects of cell coupling on beat-to-beat variability may have important implication on the electrical stability of the heart under pathological conditions characterized by altered distribution and function of gap-junctions, as previously reported for diabetes [START_REF] Nygren | Propagation of the cardiac impulse in the diabetic rat heart: reduced conduction reserve[END_REF]. Thus, altered behavior at the cellular level may have profound consequences on the entire myocardium under conditions affecting cell-to-cell coupling.

Methods presented here were applied to APs with a strictly decreasing repolarization phase, a typical feature of the electrical recovery of rodent cardiomyocytes. However, cells obtained from large mammals have APs with a prominent plateau phase [START_REF] Rota | Patch-clamp analysis in canine cardiac Purkinje cells of a novel sodium component in the pacemaker range[END_REF][START_REF] Signore | Inositol 1,4,5-trisphosphate receptors and human left ventricular myocytes[END_REF][START_REF] Sorrentino | Myocyte repolarization modulates myocardial function in aging dogs[END_REF] with repolarization not strictly decreasing. It is conceivable that segmentation of the AP repolarization represents a valid strategy for the application of our algorithms. To test this possibility, further studies are granted.

For similar reasons, methods presented here cannot be applied to AP presenting early afterdepolarizations. In a previous report [START_REF] Boutjdir | Early afterdepolarization formation in cardiac myocyte[END_REF], aberrant episodes were compared to normal APs by phase-plane plots, which exhibited oscillatory behavior in correspondence to their onset. However, the visual approach cannot be employed to quantitatively compute AP properties and does not allow the analysis of APD dynamics and statistical comparison between cell populations.

We previously documented that increasing the rate of stimulation from 2 to 6 Hz progressively exacerbated differences between the AP of control and STZ myocytes [START_REF] Meo | Reduction in Kv current enhances the temporal dispersion of the action potential in diabetic myocytes: insights from a novel repolarization algorithm[END_REF]. In the current investigation, however, focus was placed on methodological aspects of our algorithms and, therefore, the analysis of AP properties in healthy and diabetic myocytes was restricted to a single frequency of stimulation. Based on presented results, our novel methodologies represent the ideal tools to evaluate the effects of stimulation rate on the temporal evolution of the AP in normal and diseased cells, aspects that will be addressed in future investigations.

The application of our analytical approach to electrical signal obtained by optical mapping [START_REF] Beauchamp | Electrical coupling and propagation in engineered ventricular myocardium with heterogeneous expression of connexin43[END_REF][START_REF] Herron | Optical imaging of voltage and calcium in cardiac cells & tissues[END_REF] or multielectrode array methodologies [START_REF] Rossi | Ventricular activation is impaired in aged rat hearts[END_REF][START_REF] Camelliti | Adult human heart slices are a multicellular system suitable for electrophysiological and pharmacological studies[END_REF][START_REF] Kienast | Alterations of field potentials in isotropic cardiomyocyte cell layers induced by multiple endogenous pacemakers under normal and hypothermal conditions[END_REF] has the potential to allow the combination of information on repolarization properties and pattern of electrical propagation and conduction velocity. Thus, the pro-cessing of field potentials and signals obtained from voltage-sensitive molecules may lead to new modalities of analysis for ex vivo and in vitro data.

In conclusion, we have implanted methodologies allowing detailed, automatic analysis of the AP repolarization and identification of critical alterations occurring in behavior of myocytes under pathological conditions.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Upper panel: Representative APs recorded in cardiomyocytes isolated from CTRL (blue traces) and STZ (red traces) mice. Lower panel: Representation of superimposed inverted AP signals following computation of the inverse function for a train of APs obtained from a STZ myocyte. AP duration (y axis) is presented as a function of membrane voltage, expressed in percentage of repolarization (x axis).

e r 2 , under the constraint v 2 2 = 1 .

 221 Symbol '•' refers to the Hadamard elementwise product and 1I is a vector with all elements equal to 1. As formulated in Eq. (7), the problem allows the estimation of p r and v by alternated least squares minimization of the criterion function M r=1 This condition guarantees uniqueness of the solution and scalinginvariance. Based on these points, AP observations can be indicated as y r = p r • v, where p r is the polynomial p r normalized by its mean value p = 1 M r p r , and v = p • v, such that 1 M r y r = v (see Appendix). The normalization of p r by its mean value prevents the algorithm from converging to -p r and -v while obtaining the same result for the Hadamard product. Additional details related to the ESVD solution computation are provided in the Appendix.

Fig. 3 :

 3 Fig. 3: A and B. Bivariate plot of AP duration at discrete repolarization levels computed with conventional (Standard) and newly developed (Automatic, i.e. computation of the inverse function for AP signals) approaches. Parameters for the initial (START) and last AP (END) of the series for CTRL (A, n=41) and STZ (B, (n=76) myocytes are reported as mean±SEM and are fitted with linear regression. C and D. Comparison of APD properties of the initial and final AP computed with conventional (Standard) and newly developed (Automatic) approaches, for CTRL and STZ cells reported in A and B. Data is shown as median and interquartile ranges. * , p<0.05.

Fig. 4 :

 4 Fig.4: Examples of approximations of the AP profile from a STZ mouse cardiomyocyte estimated through the direct approach at fixed levels of repolarization (upper traces: 95% of repolarization; lower traces: 55% of repolarization). Blue lines, dashed red lines and solid red lines correspond to the raw data, linear regression and 5th-order polynomial approximation, respectively.

Fig. 5 :

 5 Fig. 5: Upper panel: Median value of the signed rank test for the APT DIR index computed throughout the entire repolarization using raw data (dotted lines), linear regression (dashed lines), and 5th-order polynomial approximation (solid lines) for APs obtained in CTRL (blue traces) and STZ (red traces) myocytes. Lower panel: p value of the corresponding signed rank test applied to each parameter. Black horizontal line is in correspondence of p=0.05.

Fig. 6 :

 6 Fig. 6: Statistical comparison of APT DIR obtained for CTRL and STZ groups throughout the entire repolarization range. p value of the unpaired Wilcoxon rank-sum test for APT DIR parameters computed using raw data (dotted lines), linear regression (dashed lines), and 5thorder polynomial approximation (solid lines) for the entire repolarization are reported. Black horizontal line is in correspondence of p=0.05.

Fig. 7 :Fig. 8 :Fig. 9 :

 789 Fig. 7: Upper panel: Representative example of approximation of the AP profile for a STZ myocyte during the course of stimulation (stimulation index) through dimensionless multiplicative polynomials pr calculated by ESVD method. a.u.: arbitrary units. Lower panel: Plot illustrating first-order derivatives.

Fig. 10 :Fig. 11 :Fig. 12 : 4 . Discussion and conclusions 4 . 1 .

 101112441 Fig. 10: Upper panel: Mean profile of MAPs obtained from CTRL (n=7, blue line) and STZ hearts (n=5, red line). Lower panel: Statistical comparison of MAP profile for CTRL and STZ hearts. Black horizontal line is in correspondence of p=0.05.

Acknowledgments

The authors have no disclosure to declare.

Grants

This work was supported in part by the National Institute of Health grant R01-AG-055407 and intramural resources at New York Medical College (NYMC) including funds from the NYMC Translation Science Institute. Dr Meo is funded by a Lefoulon-Delalande Foundation fellowship administered by the Institute of France.

Appendix: Computation of the ESVD solution

This section illustrates the main steps of the ESVD optimization algorithm to solve the problem: y r = p r • v + e r = B • (v1I T )θ θ θ r + e r = M v θ θ θ r + e r (A.1)

It is worth noticing that in standard SVD functions p r are replaced by scalar values a r . That would imply that if vectors y r were shifted in time, vectors v would be shifted accordingly, but SVD coefficients a r would be unaffected.

In contrast, in Eq. (A.1) functions p r and v can both capture AP temporal structure, thus allowing for a simpler data interpretation than the one obtained by conventional SVD with higher decomposition order.

Assuming that the noise samples follow a normal centered law and are uncorrelated, then the criterion to be minimized for the estimation of the vectors (A.6)

The functional J can be rewritten as:

and can be simplified to:

where diag(.) stands for the diagonal matrix. Assuming stationary conditions with respect to v, i.e., ∂J ∂v = 0, gives:

The vector v fulfilling this condition can be estimated as:

By using in sequence Eq. (A.2), Eq. (A.6) and Eq. (A.10), an alternated least square algorithm can be applied to determine the model unknowns. The initialization could be performed by choosing initial v in Eq. (A.6) as the standard mean of all the observations y r . Nevertheless, according to this formulation, the solution obtained is not unique, as multiplying v by a scalar is equivalent to multiplying p r by its inverse, yielding the same y r . This issue can be fixed by imposing the unit norm constraint over v. Since the optimization process is iterative, this normalization is directly included in the computation of Eq. (A.10).

Therefore, after convergence of the algorithm we obtain the approximation:

In order to obtain a global overview of data structure, we can compute polynomials' barycenter ṽ(i) as the center of gravity of the input observations. To this end, the expression in Eq. (A.11) can be rewritten as:

Accordingly, we obtain the necessary equality M r=1 pr (i) = α = constant. The mean value of the observations can be computed by averaging both sides of the expression in Eq. (A.12) and using the aforementioned equality:

It turns out that ṽ(i) can be computed as:

The unit norm constraint applied to ṽ enables the computation of the α value:

The computation of the modified polynomials used as descriptive features of the AP profile is finally given by the solution of the equation: p r (i)v(i) = pr (i)ṽ(i) = pr (i) 1 α r p r (i)v(i)∀i (A.16) which is:

p r (i)) (A.17)


Denys Dutykh
email: denys.dutykh@univ-savoie.fr

Visco-potential free-surface flows and long wave modelling

Keywords: free-surface flows, viscous damping, long wave models, Boussinesq equations, dissipative Korteweg-de Vries equation, bottom boundary layer PACS: 47.35.Bb, 47.35.Fg, 47.10.ad

 we presented a novel visco-potential free surface flows formulation. The governing equations contain local and nonlocal dissipative terms. From physical point of view, local dissipation terms come from molecular viscosity but in practical computations, rather eddy viscosity should be used. On the other hand, nonlocal dissipative term represents a correction due to the presence of a bottom boundary layer. Using the standard procedure of Boussinesq equations derivation, we come to nonlocal long wave equations. In this article we analyse dispersion relation properties of proposed models. The effect of nonlocal term on solitary and linear progressive waves attenuation is investigated. Finally, we present some computations with viscous Boussinesq equations solved by a Fourier type spectral method.

Introduction

Even though the irrotational theory of free-surface flows can predict successfully many observed wave phenomena, viscous effects cannot be neglected under certain circumstances. Indeed the question of dissipation in potential flows of fluid with a free surface is an important one. As stated by [START_REF] Longuet-Higgins | Theory of weakly damped stokes waves: a new formulation its physical interpretation[END_REF], it would be convenient to have equations and boundary conditions of comparable simplicity as for undamped free-surface flows. The peculiarity here lies in the fact that the viscous term in the Navier-Stokes (NS) equations is identically equal to zero for a velocity deriving from a potential. There is also a problem with boundary conditions. It is well known that the usual non-slip condition on the solid boundaries does not allow to simulate free surface flows in a finite container. Hence, some further modifications are required to permit the free surface particles to slide along the solid boundary. These difficulties were overcome in our recent work [START_REF] Dutykh | Viscous potential free-surface flows in a fluid layer of finite depth[END_REF] and [START_REF] Dias | Theory of weakly damped free-surface flows: a new formulation based on potential flow solutions[END_REF].

The effects of viscosity on gravity waves have been addressed since the end of the nineteenth century in the context of the linearized Navier-Stokes (NS) equations. It is well-known that Lamb [START_REF] Lamb | Hydrodynamics[END_REF] studied this question in the case of oscillatory waves on deep water. What is less known is that Boussinesq studied this effect as well [START_REF] Boussinesq | Lois de l'extinction de la houle en haute mer[END_REF]. In this particular case they both showed that dα dt = -2νk 2 α,

where α denotes the wave amplitude, ν the kinematic viscosity of the fluid and k = 2π/ℓ the wavenumber of the decaying wave. Here ℓ stands for the wavelength. This equation leads to the classical law for viscous decay, namely

α(t) = α 0 e -2νk 2 t . (2)
Let us consider a simple numerical application with g = 9.8 m/s 2 , ℓ = 3 m and molecular viscosity ν = 10 -6 m 2 /s. According to the formula (2), this wave will take t 0 ≈ 8 × 10 4 s or about one day before losing one half of its amplitude. This wave will attain velocity

c p = g k ≈ 2.16 m/s
and travel the distance equal to L = c p • t 0 ≈ 170 km. This estimation is exaggerated since the classical result of Boussinesq and Lamb does not take into account energy dissipation in the bottom boundary layer. We will discuss the question of linear progressive waves attenuation in Section 5. Another point is that the molecular viscosity ν should be replaced by eddy viscosity ν t which is more appropriate in most practical situations (see Remark 6 for more details).

The importance of viscous effects for water waves has been observed in various experimental studies. For example, in [START_REF] Zabusky | Shallow water waves, the Korteweg-de Vries equation and solitons[END_REF] one can read . . . However, the amplitude disagrees somewhat, and we suppose that this might be due to the viscous dissipation. . . [START_REF] Wu | Long waves in ocean and coastal waters[END_REF] also mentions this drawback of the classical water wave theory:

. . . the peak amplitudes observed in the experiments are slightly smaller than those predicted by the theory. This discrepancy can be ascribed to the neglect of the viscous effects in the theory. . .

Another example is the conclusion of [START_REF] Bona | An evaluation of a model equation for water waves[END_REF]:

. . . it was found that the inclusion of a dissipative term was much more important than the inclusion of the nonlinear term, although the inclusion of the nonlinear term was undoubtedly beneficial in describing the observations. . .

Another source of dissipation is due to bottom friction. An accurate computation of the bottom shear stress τ is crucial for calculating sediment transport fluxes. Consequently, the predicted morphological changes will greatly depend on the chosen shear stress model. Traditionally, this quantity is modelled by a Chézy-type law

τ | z=-h = C f ρ| u h | u h ,
where u h = u(x, y, -h, t) is the fluid velocity at the bottom, C f is the friction coefficient and ρ is the fluid density. Two other often used laws can be found in [DD07a, Section 4.3], for example. One problem with this model is that τ and u h are in phase. It is well known [LSVO06] that in the case of a laminar boundary layer, the bottom stress τ | z=-h is π 4 out of phase with respect to the bottom velocity.

Water wave energy can be dissipated by different physical mechanisms. The research community agrees at least on one point: the molecular viscosity is unimportant. Now let us discuss more debatable statements. For example if we take a tsunami wave and estimate its Reynolds number, we find Re ≈ 10 6 . So, the flow is clearly turbulent and in practice it can be modelled by various eddy viscosity models. On the other hand, in laboratory experiments the Reynolds number is much more moderate and sometimes we can neglect this effect. When nonbreaking waves feel the bottom, the most efficient mechanism of energy dissipation is the bottom boundary layer. This is the focus of our paper. We briefly discuss the free surface boundary layer and explain why we do not take it into account in this study. Finally, the most important (and the most challenging) mechanism of energy dissipation is wave breaking. This process is extremely difficult from the mathematical but also the physical and numerical points of view since we have to deal with multivalued functions, topological changes in the flow and complex turbulent mixing processes. Nowadays the practitioners can only be happy to model roughly this process by adding adhoc dissipative terms when the wave becomes steep enough.

In this work we keep the features of undamped free-surface flows while adding dissipative effects. The classical theory of viscous potential flows is based on pressure and boundary conditions corrections [START_REF] Joseph | The dissipation approximation and viscous potential flow[END_REF] due to the presence of viscous stresses. We present here another approach.

Currently, potential flows with ad-hoc dissipative terms are used for example in direct numerical simulations of weak turbulence of gravity waves [DKZ03, DKZ04, ZKPD05]. There have also been several attempts to introduce dissipative effects into long wave modelling [START_REF] Mei | The applied dynamics of ocean surface waves[END_REF][START_REF] Dutykh | Dissipative Boussinesq equations[END_REF][START_REF] Chen | Long-time asymptotic behaviour of dissipative Boussinesq systems[END_REF]Kha97]. We would like to underline that the last paper [Kha97] contains also a nonlocal dissipative term in time.

The present work is a direct continuation of the recent studies [START_REF] Dias | Theory of weakly damped free-surface flows: a new formulation based on potential flow solutions[END_REF][START_REF] Dutykh | Viscous potential free-surface flows in a fluid layer of finite depth[END_REF]. In [START_REF] Dias | Theory of weakly damped free-surface flows: a new formulation based on potential flow solutions[END_REF] the authors considered periodic waves in infinite depth and derivation was done in two-dimensional (2D) case, while in [START_REF] Dutykh | Viscous potential free-surface flows in a fluid layer of finite depth[END_REF] we removed these two hypotheses and all the computations are done in 3D. This point is important since the vorticity structure is more complicated in 3D. In other words we considered a general wavetrain on the free surface of a fluid layer of finite depth. As a result we obtained a formulation which contains a nonlocal term in the bottom kinematic condition. Later we discovered that the nonlocal term in exactly the same form was derived in [LO04]. The inclusion of this term is natural since it represents the correction to potential flow due to the presence of a boundary layer. Moreover, this term is predominant since its magnitude scales with O( √ ν), while other terms in the free-surface boundary conditions are of order O(ν). The importance of this effect was pointed out in the classical literature on the subject [Lig78]:

. . . Bottom friction is the most important wherever the water depth is substantially less than a wavelength so that the waves induce significant horizontal motions near the bottom; the associated energy dissipation takes place in a boundary layer between them and the solid bottom. . . This quotation means that this type of phenomenon is particularly important for shallow water waves like tsunamis, for example [START_REF] Dias | Extreme Man-Made and Natural Hazards in Dynamics of Structures, chapter Dynamics of tsunami waves[END_REF]Dut07]. Here we present several numerical computations based on the newly derived governing equations and analyse dispersion relation properties.

We would like to mention here a paper of N. Sugimoto [START_REF] Sugimoto | Burgers equation with a fractional derivative; hereditary effects on nonlinear acoustic waves[END_REF]. The author considered initial-value problems for the Burgers equation with the inclusion of a hereditary integral known as the fractional derivative of order 1 2 . The form of this term was justified in previous works [Sug89, Sug90]. Note, that from fractional calculus point of view our nonlocal term (6) is also a half-order integral.

Other researchers have obtained nonlocal corrections but they differ from ours [KM75]. This discrepancy can be explained by a different scaling chosen by Kakutani & Matsuuchi in the boundary layer. Consequently, their governing equations contain a nonlocal term in space. The performance of the present nonlocal term (6) was studied in [LSVO06]. The authors carried out in a wave tank a set of experiments, analyzing the damping and shoaling of solitary waves. It is shown that the viscous damping due to the bottom boundary layer is well represented by the nonlocal term. Their numerical results fit well with the experiments. The model not only properly predicts the wave height at a given point but also provides a good representation of the changes on the shape and celerity of the soliton. We can conclude that the experimental study by P. Liu et al. [LSVO06] validates this theory.

The present article is organized as follows. In Section 2 we estimate the rate of viscous dissipation in different regions of the fluid domain. Then, we present basic ideas of derivation and come up with visco-potential free-surface flows formulation. At the end of Section 3 we give corresponding long wave models: nonlocal Boussinesq and KdV equations. Section 4 is completely devoted to the analysis of linear dispersion relation of complete and long wave models introduced in previous section. Last two sections deal with linear progressive and solitary waves attenuation respectively. Finally, the paper is ended by some conclusions and perspectives.

Anatomy of dissipation

In this section we briefly discuss the contribution of different flow regions into water wave energy dissipation. We conventionally [START_REF] Mei | The applied dynamics of ocean surface waves[END_REF] divide the flow into three regions illustrated on Figure 1. On this figure S f and S b stand for free surface and bottom respectively. Then, R i , R f and R b denote the interior region, free surface and bottom boundary layers. In order to make some estimates we introduce the notation which will be used in this section: µ is the dynamic viscosity, δ = O( √ µ) is the boundary layer thickness, t 0 is the characteristic time, a 0 is the characteristic wave amplitude and ℓ is the wavelength.

O x z S f S b R f R b R i
We assume that the flow is governed by the incompressible Navier-Stokes equations:

∇ • u = 0, ∂ u ∂t + u • ∇ u + 1 ρ ∇p = g + 1 ρ ∇ • τ ,
where τ is the viscous stress tensor

τ ij = 2µε ij , ε ij = 1 2 ∂u i ∂x j + ∂u j ∂x i .
We multiply the second equation by u and integrate over the domain Ω with boundary ∂Ω to get the following energy balance equation:

1 2 Ω ∂ ∂t ρ| u| 2 dΩ + 1 2 ∂Ω ρ| u| 2 u • n dσ = = ∂Ω -pI + τ n • u dσ + Ω ρ g • u dΩ - 1 2µ Ω τ : τ dΩ T .
In this identity each term has a precise physical meaning. The left-hand side is the total rate of energy change in Ω. The second term is the flux of energy across the boundary. On the right-hand side, the first integral represents the rate of work by surface stresses acting on the boundary. The second integral is the rate of work done by the gravity force throughout the volume, and the third integral T is the rate of viscous dissipation. We focus our attention on the last term T . We estimate the order of magnitude of the rate of dissipation in various regions of the fluid.

We start by the interior region R i . Outside the boundary layers, it is reasonable to expect that the rate of strain is dominated by the irrotational part of the velocity whose scale is a 0 t 0 and the length scale is the wavelength ℓ. The energy dissipation rate is then

O T R i ∼ 1 µ µ a t 0 ℓ 2 • ℓ 3 = µ a t 0
of u dominates the strain rate, so that

O T R b ∼ 1 µ µ a t 0 δ 2 • δℓ 2 = µ δ aℓ t 0 2 ∼ O(µ 1 2 ) .
A free surface boundary layer also exists. Its importance depends on the free surface conditions. Consider first the classical case of a clean surface. The stress is mainly controlled by the potential velocity field which is of the same order as in the main body of the fluid. Because of the small volume O(δℓ 2 ) the rate of dissipation in the free surface boundary layer is only

O T R f ∼ 1 µ µ a t 0 ℓ 2 • δℓ 2 = µδ a t 0 2 ∼ O(µ 3 2 ) .
From the physical point of view it is weaker, since only the zero shear stress condition on the free surface is required.

Another extreme case is when the free surface is heavily contaminated, for example, by oil slicks. The stress in the free surface boundary layer can then be as great as in the boundary layer near a solid wall. In the present study we do not treat such extreme situations and the surface contamination is assumed to be absent.

The previous scalings suggest the following diagram which represents the hierarchy of dissipative terms:

O µ 1 2 R b ֒→ O(µ) R i ֒→ O µ 3 2 R f ֒→ . . .
It is clear that the largest energy dissipation takes place inside the bottom boundary layer. We take into account only two first phenomena from this diagram. Consequently, all dissipative terms of order O(µ

2 ) and higher will be neglected.

Remark 1 In laboratory experiments, surface waves are confined by side walls as well. The rate of attenuation due to the side-wall boundary layers was computed in [START_REF] Mei | The damping of surface gravity waves in a bounded liquid[END_REF]. For simplicity, in the present study we consider an unbounded in horizontal coordinates domain (see Figure 1).

Derivation

Consider the linearized 3D incompressible NS equations describing free-surface flows in a fluid layer of uniform depth h:

∂ v ∂t = - 1 ρ ∇p + ν∆ v + g, ∇ • v = 0, (3)
with v the velocity vector, p the pressure, ρ the fluid density and g the acceleration due to gravity. We represent v = (u, v, w) in the form of the Helmholtz-Leray decomposition:

v = ∇φ + ∇ × ψ, ψ = (ψ 1 , ψ 2 , ψ 3 ). ( 4
)
After substitution of the decomposition (4) into (3), one notices that the equations are verified provided that the functions φ and ψ satisfy the following equations:

∆φ = 0, φ t + p -p 0 ρ + gz = 0, ∂ ψ ∂t = ν∆ ψ.
Next we discuss the boundary conditions. We assume that the velocity field satisfies the conventional no-slip condition at the bottom v| z=-h = 0, while at the free surface we have the usual kinematic condition, which can be stated as η t + v • ∇η = w. After linearization it becomes simply η t = w.

Dynamic condition states that the forces must be equal on both sides of the free surface:

[ σ • n] ≡ -(p -p 0 ) n + τ • n = 0 at z = η(x, t) ,
where σ is the stress tensor, [f ] denotes the jump of a function f across the free surface, n is the normal to the free surface and τ the viscous part of the stress tensor σ.

The basic idea consists in expressing the vortical part of the velocity field ∇ × ψ in terms of the velocity potential φ and the free surface elevation η using differential or pseudodifferential operators. In this section we just show final results while the details of computation can be found in [START_REF] Dutykh | Viscous potential free-surface flows in a fluid layer of finite depth[END_REF]Dut07].

Let us begin by the free-surface kinematic condition

∂η ∂t = w ≡ ∂φ ∂z + ∂ψ 2 ∂x - ∂ψ 1 ∂y , z = 0.
Using the absence of tangential stresses on the free surface, one can replace the rotational part in the kinematic boundary condition:

η t = φ z + 2ν∆η.
In order to account for the presence of viscous stresses, we have to modify the dynamic free-surface condition as well. This is done using the balance of normal stresses at the free surface:

σ zz = 0 at z = 0 ⇒ p -p 0 = 2ρν ∂w ∂z ≡ 2ρν ∂ 2 φ ∂z 2 + ∂ 2 ψ 2 ∂x∂z - ∂ 2 ψ 1 ∂y∂z .
In [START_REF] Dutykh | Viscous potential free-surface flows in a fluid layer of finite depth[END_REF] it is shown that

∂ 2 ψ 2 ∂x∂z -∂ 2 ψ 1 ∂y∂z = O(ν 1
2 ), so Bernoulli's equation becomes

φ t + gη + 2νφ zz + O(ν 3
2 ) = 0.

(5) Since we only consider weak dissipation (ν ∼ 10 -6 -10 -3 m 2 /s), we neglect terms of order o(ν).

The second step in our derivation consists in introducing a boundary layer correction at the bottom. In the present study we assume the boundary layer to be laminar. The question of turbulent boundary layer was investigated in [Liu06]. Hence, the bottom boundary condition becomes

∂φ ∂z z=-h = ν π t 0 ∇ 2 x φ 0 | z=-h √ t -τ dτ = - ν π t 0 φ 0zz | z=-h √ t -τ dτ. (6)
One recognizes on the right-hand side a half-order integral operator. The last equation shows that the effect of the diffusion process in the boundary layer is not instantaneous. The result is cumulative and it is weighted by (t -τ ) -1 2 in favour of the present time.

Summarizing the developments made above and generalizing our equations by including nonlinear terms, we obtain a set of viscous potential free-surface flow equations:

∆φ = 0, ( x, z) ∈ Ω = R 2 × [-h, η] (7) η t + ∇η • ∇φ = φ z + 2ν∆η, z = η (8) φ t + 1 2 |∇φ| 2 + gη = -2νφ zz , z = η (9) φ z = - ν π t 0 φ zz √ t -τ dτ, z = -h. ( 10
)
At the present stage, the addition of nonlinear terms is rather a conjecture. However, a recent study by Liu et al. [START_REF] Liu | Boundary layer flow and bed shear stress under a solitary wave[END_REF] suggests that this conjecture is rather true. The authors investigated the importance of nonlinearity in the bottom boundary layer for a solitary wave solution. They came to the conclusion that "the nonlinear effects are not very significant".

Using this weakly damped potential flow formulation and described in previous works procedure [DD07a, Section 4] of Boussinesq equations derivation, one can derive the following system of equations with horizontal velocity u θ defined at the depth z θ = -θh, 0 ≤ θ ≤ 1:

η t + ∇ • ((h + η) u θ ) + h 3 θ 2 2 -θ + 1 3 ∇ 2 (∇ • u θ ) = 2ν∆η + ν π t 0 ∇ • u θ √ t -τ dτ, (11)
u θt + 1 2 ∇| u θ | 2 + g∇η -h 2 θ 1 - θ 2 ∇(∇ • u θt ) = 2ν∆ u θ . ( 12
)
For simplicity, in this study we present governing equations only on the flat bottom, but generalization can be done for general bathymetry. Khabakhpashev [Kha87], Liu & Orfila [LO04] also derived a similar set of Boussinesq equations in terms of depth-averaged velocity. Our equations ( 11) -( 12) have local dissipative terms in addition and they are formulated in terms of the velocity variable defined at arbitrary water level that is beneficial for dispersion relation properties.

Total energy decay in visco-potential flow

Let us consider a fluid layer infinite in horizontal coordinates x = (x, y), bounded below by the flat bottom z = -h and above by the free surface. Total energy of water waves is given by the following formula:

E = R 2 η -h 1 2 |∇φ| 2 dzd x + g 2 R 2 η 2 d x.
We would like to mention here that Zakharov showed [Zak68] this expression to be the Hamiltonian for classical water wave problem with suitable choice of canonical variables: η and ψ := φ( x, z = η, t).

In this section we are interested in the evolution of the total energy E with time. This question is investigated by computing the derivative dE dt with respect to time t. Obviously, when one considers the classical potential free surface flow formulation [START_REF] Lamb | Hydrodynamics[END_REF], we have dE dt ≡ 0, since no dissipation is introduced into the governing equations. We performed the computation of dE dt in the framework of the visco-potential formulation and give here the final result:

dE dt = ν π R 2 φ t | z=-h t 0 φ zz | z=-h √ t -τ dτ d x + 2ν R 2 ∂ t ( ∇φ| z=η • ∇η) d x.
In this identity, the first term on the right hand side comes from the boundary layer and is predominant in the energy decay since its magnitude scales with O( √ ν). The second term has its origins in free surface boundary conditions.

Its magnitude is O(ν), thus it has less important impact on the energy balance. This topic will be investigated further in future studies.

Dissipative KdV equation

In this section we derive a viscous Korteweg-de Vries (KdV) equation from just obtained Boussinesq equations ( 11), (12). Since KdV-type equations model only unidirectional wave propagation, our attention is naturally restricted to 1D case. In order to perform asymptotic computations, all the equations have to be switched to nondimensional variables as it is explained in [DD07a, Section 2]. We find the velocity variable u in this form:

u θ = η + εP + µ 2 Q + . . .
where ε, µ are nonlinearity and dispersion parameters respectively (see [START_REF] Dutykh | Dissipative Boussinesq equations[END_REF] for their definition), P and Q are unknown at the present moment. Using the methods similar to those used in [DD07a, Section 6.1], one can easily show that

P = - 1 4 η 2 , Q = θ - 1 6 - θ 2 2 η xx .
This result immediately yields the following asymptotic representation of the velocity field

u θ = η - 1 4 εη 2 + µ 2 θ - 1 6 - θ 2 2 η xx + . . . (13)
Substituting the last formula (13) into equation ( 11) and switching again to dimensional variables, one obtains this viscous KdV-type equation: . . . it may be concluded that our modified K-dV equation can describe the observed wave behaviours except the fact that the phase shift obtained by the calculations is not confirmed by their experiments.

η t + g h (h + 3 2 η)η x + 1 6 h 3 η xxx - ν π t 0 η x √ t -τ dτ = 2νη xx . (14
Excellent performance of the model ( 14) with respect to experiments was shown in [LSVO06].

Dispersion relation of complete and Boussinesq nonlocal equations

Interesting information about the governing equations can be obtained from the linear dispersion relation analysis. In this section we are going to analyse the new set of equations ( 7)-(10) for the complete water wave problem and the corresponding long wave asymptotic limit (11), (12).

To simplify the computations, we consider the two-dimensional problem. The generalization to higher dimensions is straightforward and is performed by replacing the wavenumber k by its modulus | k| in vectorial case. Traditionally the governing equations are linearized and the bottom is assumed to be flat. The last hypothesis is made throughout this study. After all these simplifications the new set of equations becomes

φ xx + φ zz = 0, (x, z) ∈ R × [-h, 0], ( 15
) η t = φ z + 2νη xx , z = 0, ( 16
) φ t + gη + 2νφ zz = 0, z = 0, (17)
φ z + ν π t 0 φ zz √ t -τ dτ = 0, z = -h. ( 18
)
The next classical step consists in finding solutions of the special form

φ(x, z, t) = ϕ(z)e i(kx-ωt) , η(x, t) = η 0 e i(kx-ωt) . ( 19
)
From continuity equation (15) we can determine the structure of the function ϕ(z):

ϕ(z) = C 1 e kz + C 2 e -kz .
Altogether we have three unknown constants2 C = (C 1 , C 2 , η 0 ) and three boundary conditions ( 16)-( 18) which can be viewed as a linear system with respect to C:

M C = 0. ( 20
)
The matrix M has the following elements

M =        k -k iω -2νk 2 iω -2νk 2 iω -2νk 2 -g e -kh 1 + kF (t, ω) e kh -1 + kF (t, ω) 0       
where the function F (t, ω) is defined in the following way:

F (t, ω) := ν π t 0 e iω(τ )(t-τ ) √ t -τ dτ. (21)
In order to have nontrivial solutions of ( 15)-( 18), the determinant of the system (20) has to be equal to zero det M = 0. It gives us a relation between ω and wavenumber k. This relation is called the linear dispersion relation:

D(ω, k) := (iω -2νk 2 ) 2 + gk tanh(kh) -kF (t, ω) (iω -2νk 2 ) 2 tanh(kh) + gk ≡ 0. ( 22
)
A similar procedure can be followed for Boussinesq equations (11), (12). We do not give here the details of the computations but only the final result:

D b (ω, k) := (iω -2νk 2 ) 2 + b(kh) 2 iω(iω -2νk 2 ) + ghk 2 (1 -a(kh) 2 ) -gk 2 F (t, ω) ≡ 0, ( 23
)
where we introduced the following notation:

a := θ 2 2 -θ + 1 3 , b := θ 1 -θ 2 .
Unfortunately, the relations D(ω, k) ≡ 0 and D b (ω, k) ≡ 0 cannot be solved analytically to give an explicit dependence of ω on k as for the classical water wave problem. Consequently, we apply a quadrature formula to discretize the nonlocal term F (t n , ω), where t n denotes the discrete time variable. The resulting algebraic equation with respect to ω(t n ; k) is solved analytically.

Remark 2 Contrary to the classical water wave problem and, by consequence, standard Boussinesq equations (their dispersion relation can be found in [DD07a, Section 3.2], for example) where the dispersion relation does not depend on time ω 2 -gk tanh(kh) ≡ 0, (24)

here we have additionally the dependence of ω(k; t) on time t as a parameter. It is a consequence of the presence of the nonlocal term in time in the bottom boundary condition (10). Physically it means that the boundary layer "remembers" the flow history.

Remark 3 There is one subtle point in the derivation presented above. In fact, all computations were performed as if the frequency ω were independent from time. Our final result shows that time t appears explicitly in the dispersion relations ( 22), (23). Developments made above make sense under the assumption of slow variation of ω with time t. This statement can be written in mathematical form ∂ω ∂t ≪ 1. It is rather a conjecture here and will be examined in future studies. We had to make this assumption in order to avoid complicated integro-differential equations and, consequently, simplify the analysis.

Analytical limit for infinite time

In the previous section we showed that the dispersion relation of our viscopotential formulation is time dependent. It is natural to ask what happens when time evolves. Here we compute the limiting state of the dispersion curves (22), (23) as t → +∞. Namely, we will take this limit in equations ( 22), (23) assuming, of course, that it exists

∃ ω ∞ (k) := lim t→+∞ ω t (k).
Time t comes in dispersion relations through the argument of the function F (t, ω) defined in (21). Its limit can be easily computed to give

lim t→+∞ F (t, ω) = ν πω ∞ +∞ 0 e ip √ p dp = ν ω ∞ e i π 4 .
Now we are ready to write down the final results:

D(ω ∞ , k) := (iω ∞ -2νk 2 ) 2 + gk tanh(kh) - ν ω ∞ ke i π 4 (iω ∞ -2νk 2 ) 2 tanh(kh) + gk ≡ 0, D b (ω ∞ , k) := (iω ∞ -2νk 2 ) 2 + b(kh) 2 iω ∞ (iω ∞ -2νk 2 ) + ghk 2 (1 -a(kh) 2 ) -gk 2 ν ω ∞ e i π 4 ≡ 0.
In order to solve numerically nonlinear equation D(ω ∞ , k) = 0 (or D b (ω ∞ , k) = 0 when one is interested in Boussinesq equations) with respect to ω ∞ , we apply a Newton-type method. The iterations converge very quickly since we use analytical expressions for the Jacobian dD dω∞ ( dD b dω∞ , correspondingly). Derivative computation is straightforward. Limiting dispersion curves are plotted (see Figure 7) and discussed in the next section.

Discussion

Numerical snapshots of the nonclassical dispersion relation3 at different times for complete and Boussinesq equations are given on Figures (2)-( 6). The value of the eddy viscosity ν is taken from Table 1 and we consider a one meter depth fluid layer (h = 1 m). We will try to make several comments on the obtained results.

Remark 4 Recall that these snapshots were obtained under the assumption that ω is slowly varying in time. The validity of this approximation is examined and discussed in [START_REF] Dutykh | Group and phase velocities in the free-surface viscopotential flow: new kind of boundary layer induced instability[END_REF].

Just at the beginning (when t = 0), there is no effect of the nonlocal term. This is why on Figure 2 new and classical4 curves are superimposed. With no surprise, the phase velocity of Boussinesq equations represents well only long waves limit (let us say up to kh ≈ 2). When time evolves, we can see that the main effect of nonlocal term consists in slowing down long waves (see Figures 345). Namely, in the vicinity of kh = 0 the real part of the phase velocity is slightly smaller with respect to the classical formulation. From physical point of view this situation is comprehensible since only long waves "feel" the bottom and, by consequence, are affected by bottom boundary layer. On the other hand, the imaginary part of the phase velocity is responsible for the wave amplitude attenuation. The minimum of Im c p (k) in the region of long waves indicates that there is a "preferred" wavelength which is attenuated the most. In the range of short waves the imaginary part is monotonically decreasing. In practice it means that high-frequency components are damped by the model. This property can be advantageous in numerics, for example. On Figure 6 we depicted the real part of c p (k) with zoom made on long and moderate waves. The reader can see that nonlocal full and Boussinesq equations have similar behaviour in the vicinity of kh = 0. Now let us discuss the limiting state of phase velocity curves as t → +∞. It is depicted on Figure 7. One can see singular behaviour in the vicinity of zero. This situation is completely normal since very long waves are highly affected by bottom boundary layer.

We would like to comment more on the behaviour of curves on Figure 7 since it is not easy to distinguish them with graphical resolution. We will concentrate on the upper image because everything is clear below with imaginary part. In the vicinity of kh = 0 we have the superposition of nonlocal models (complete set of equations ( 7) -(10) and Boussinesq equations ( 11) -( 12)). When we gradually move to short waves, we have the superposition of complete classical and nonlocal (7) -(10) water wave problems. Meanwhile, Boussinesq system gives slightly different phase velocity for kh ≥ 3. This is comprehensible since we cannot simplify considerably the problem and have uniformly good approximation everywhere. Various Boussinesq systems are designed to reproduce the behaviour of long waves.

Attenuation of linear progressive waves

In this Section we investigate the damping rate of linear progressive waves. Thus, the first step consists in linearizing dissipative KdV equation (14) to obtain the following nonlocal Airy equation:

η t + g h hη x + 1 6 h 3 η xxx - ν π t 0 η x √ t -τ dτ = 2νη xx (25)
In other words, we can say that we restrict our attention only to small amplitude waves. Now we make the next assumption. We look for a particular form of the solutions:

η(x, t) = A(t)e ikξ , ξ = x -ght. ( 26
)
where k is the wavenumber and A(t) is called the complex amplitude, since |η(x, t)| = |A(t)|. Integro-differential equation governing the temporal evolution of A(t) can be easily derived by substituting the special representation (26) into linearized KdV equation (25):

dA dt - i 6 g h (kh) 3 A(t) + 2νk 2 A(t) -ik gν πh t 0 A(τ ) √ t -τ dτ = 0. (27)
In our applications we are rather interested in temporal evolution of the absolute value |A(t)|. It is straightforward to derive the governing equation for the squared wave amplitude:

d|A| 2 dt + 4νk 2 |A(t)| 2 -ik gν πh t 0 Ā(t)A(τ ) -A(t) Ā(τ ) √ t -τ dτ = 0.
If we denote by A r (t) and A i (t) real and imaginary parts of A(t) respectively, the last equation can be further simplified:

d|A| 2 dt + 4νk 2 |A(t)| 2 + 2k gν πh t 0 A r (t)A i (τ ) -A i (t)A r (τ ) √ t -τ dτ = 0. ( 28
)
Just derived integro-differential equation represents a generalisation to the classical equation (1) by Boussinesq [START_REF] Boussinesq | Lois de l'extinction de la houle en haute mer[END_REF] and Lamb [START_REF] Lamb | Hydrodynamics[END_REF] for the peri-odic, linear wave amplitude evolution in a viscous fluid. We recall that novel integral term is a direct consequence of the bottom boundary layer modelling.

Unfortunately, equation (28) cannot be used directly for numerical computations since we need to know the following combination of real and imaginary parts A r (t)A i (τ ) -A i (t)A r (τ ) for τ ∈ [0, t]. It represents a new and nonclassical aspect of the present theory.

Equation (28) allows us to discuss the relative importance of local and nonlocal dissipative terms for long waves. In fact, when we consider the deep-water approximation, only local dissipative terms are present in the governing equations [START_REF] Dias | Theory of weakly damped free-surface flows: a new formulation based on potential flow solutions[END_REF]. On the other hand, in shallow waters the integral term is predominant. It means that there is an intermediate depth where both dissipative terms have equal magnitude. This depth can be estimated when one switches to dimensionless form of the equation ( 28). Comparing the coefficients in front of dissipative terms gives the following transcendental equation for the "critical" depth h * :

h * = g 4πωνk 2
where ω is the characteristic wave frequency.

In numerical computations it is advantageous to integrate exactly local terms in equation ( 27). It is done by making the following change of variables:

A(t) = e -2νk 2 t e i 6 √ g h (kh) 3 t Ã(t).
One can easily show that new function Ã(t) satisfies the following equation:

d à dt = ik gν πh t 0 e 2νk 2 (t-τ ) e -i 6 √ g h (kh) 3 (t-τ ) √ t -τ Ã(τ ) dτ
On Figure 8 we plot a solution of integro-differential equation ( 27). All parameters related to this case are given in Table 1. These values were chosen to simulate a typical tsunami in Indian Ocean [START_REF] Dias | Extreme Man-Made and Natural Hazards in Dynamics of Structures, chapter Dynamics of tsunami waves[END_REF]. We have to say that the wave amplitude damping is entirely due to the dissipation in boundary layer since local terms are unimportant for sufficiently long waves. It means that classical formula (2) gives almost constant value α 0 of the amplitude α(t) on the time scale of several hours, since the factor 2νk 2 is of order ≈ 10 -11 for parameters given in Table 1.

Solitary wave propagation

In this section we would like to show the effect of nonlocal term on the solitary wave attenuation. For simplicity, we will consider wave propagation in a 1D channel. The question of the bottom shear stress effect on the solitary wave propagation was considered for the first time in [Keu48].

For numerical integration of equations ( 11), (12) we use the same Fouriertype spectral method that was described in [DD07a, Section 5]. Obviously this method has to be slightly adapted because of the presence of nonlocal in time term. We have to say that this term necessitates the storage of ∇ • u (n) at previous time steps. Hence, long time computations can be memory consuming.

The values of all dimensionless parameters are given in Table 2. Dimensionless viscosity ν is related to other physical parameters in the following way:

ν 2 = ν ℓ √ gh ,
where ν is kinematic viscosity, ℓ is the characteristic wavelength and h is the typical depth.

Remark 5 From numerical point of view, the integral term

t 0 φzz( x,-h,τ ) √ t-τ
dτ can pose some problems in the vicinity of the upper limit τ = t. Probably the best way to deal with this problem is to separate the integral in two parts:

t 0 φ zz ( x, -h, τ ) √ t -τ dτ = t-δ 0 φ zz ( x, -h, τ ) √ t -τ dτ + t t-δ φ zz ( x, -h, τ ) √ t -τ dτ, δ > 0.
The first integral can be computed in a usual way without any special care.

Then one applies to the second integral a special class of Gauss quadrature formulas with weighting function (t -τ ) -1 2 .

But there is another well-known trick that we describe here. This technique can be implemented in simpler way. We rewrite our integral in the following way

t 0 φ zz ( x, -h, τ ) √ t -τ dτ = t 0 φ zz ( x, -h, t) √ t -τ dτ + t 0 φ zz ( x, -h, τ ) -φ zz ( x, -h, t) √ t -τ dτ.
The first integral in the right hand side can be evaluated analytically while the second one does not contain any singularity under the assumption of differentiability of τ → φ zz ( x, -h, τ ) at τ = t:

t 0 φ zz ( x, -h, τ ) √ t -τ dτ = 2 √ tφ zz ( x, -h, t) + t 0 φ zz ( x, -h, τ ) -φ zz ( x, -h, t) √ t -τ dτ.
Remark 6 What is the value of ν to be taken in numerical simulations? There is surprisingly little published information of this subject. What is clear is that the molecular diffusion is too small to model true viscous damping and one should rather consider the eddy viscosity parameter. Some interesting information on this subject can be found in [START_REF] Tuck | Sea wave pattern evaluation[END_REF]:

We have spent a considerable amount of time and effort seeking further published information on viscous effects on ship waves.

θ z θ = -θh 1 - √ 5/5
x 0 soliton center at t = 0 -1.5

Table 2 Values of the dimensionless parameters used in the numerical computations.

. . . we reiterate that a viscosity of about ν = 0.005 m 2 /s gave reasonable agreement with longitudinal cut results (including apparent damping of transverse waves).

In another famous paper [START_REF] Bona | An evaluation of a model equation for water waves[END_REF] one can find:

. . . Such a decay rate leads to a value for µ in (M * ) of 0.014. . .

In their work µ is the coefficient in front of the dissipative term in BBM equation:

η t + η x + 3 2 ηη x -µη xx - 1 6 η xxt = 0.
It is important to underline that this equation is written in dimensionless variables. Thus, the value reported in that study has to be rescaled with respect to other physical parameters.

Liu and Orfilla [LO04] report the value of eddy viscosity ν = 0.001 m 2 /s.

If we summarize all the remarks made above, we can conclude that the value of the order 10 -3 -10 -2 m 2 /s should give reasonable results.

Approximate solitary wave solution

In order to provide an initial condition for equations (11), (12), we are going to obtain an approximate solitary wave solution for nondissipative 1D version of these equations over the flat bottom: Then, we apply the same approach as in Section 3.2 or in [DD07a, Section 6.1]. We do not provide the computations here since they are simple and can be done without any difficulties. The final result is the following:

η t + (1 + εη)u x + µ 2 θ 2 2 -θ + 1 3 u xxx = 0, u t + η x + ε 2 (u 2 ) x -µ 2 θ 1 - θ 2 u xxt = 0.
η(x, t) = 2(c -1) ε sech 2 6(c -1) 2µ (x + x 0 -ct)
and the velocity is given by formula (13). In the numerical results presented here, we use η(x, 0) and u(x, 0) as initial conditions.

Discussion

On Figures 9-11 we present three curves. They depict the free surface elevation according to three different formulations. The first corresponds to classical Boussinesq equations without dissipation. The second one to dissipative system with differential or local terms (for example, ν∆ u in momentum conservation equation) and the third curve corresponds to the new set of equations (11), (12). On Figure 12 we made a zoom on the soliton crest.

It can be seen that Boussinesq equations with nonlocal term provide stronger attenuation of the amplitude. In the same time, as it was shown in the previous section, this nonlocal term slightly slows down the solitary wave.

In order to show explicitly the rate of amplitude experiments show that their amplitude decreases when N → ∞. This result shows one more time that nonlocal model provides stronger damping properties. One can have the impression that the amplitude decays linearly but it is only an impression because of (2). This behaviour for small t can be explained by simple Taylor expansion which is valid when νk 2 t ≪ 1:

α(t) = α 0 e -

Conclusion

In the present article we used a novel visco-potential formulation proposed recently in [LO04, [START_REF] Dutykh | Viscous potential free-surface flows in a fluid layer of finite depth[END_REF]Dut07]. This formulation contains the nonlocal term in the kinematic bottom boundary condition. This term should be considered as a boundary layer correction at the bottom. In modelling viscous effects this term plays the main role, since its magnitude is O( √ ν). Recently, important progress was made in numerical computation of the nonlocal term [TL07]. Thus, we can hope to have its implementation in operational numerical codes. It is interesting to note that the boundary layer effect is not instantaneous but rather cumulative. The flow history is weighted by (t -τ ) -1 2 in favour of the current time t. As pointed out in [LO04], this nonlocal term is essential to have an accurate estimation of the bottom shear stress based on the calculated wave field above the bed. Then, this information can be used for calculating sediment-bedload transport fluxes and, in turn, morphological changes.

A long wave asymptotic limits (Boussinesq and KdV equations) were derived from this new formulation. The dispersion relation was described. Due to the presence of special functions, we cannot obtain a simple analytical dependence of the frequency ω on the wavenumber k as in classical equations. Consequently the dispersion curve was obtained numerically by a Newton-type method. We made a comparison between the phase velocity of the complete visco-potential problem and the corresponding Boussinesq equations. The dispersion relation of new formulation is shown to be time dependent and this property is not classical. It comes from the memory effect of the boundary layer. We computed analytically the limiting dispersion curve ω ∞ (k) as t → +∞.

The effect of the nonlocal term on the solitary wave attenuation was investigated numerically with a Fourier-type spectral method. We showed that it provides much stronger damping (see Figure 13). An equation describing the amplitude evolution of linear progressive waves was derived (27). This result includes bottom boundary layer correction and generalizes classical formula (1) by Boussinesq [START_REF] Boussinesq | Lois de l'extinction de la houle en haute mer[END_REF] and Lamb [START_REF] Lamb | Hydrodynamics[END_REF].

The present study opens new exciting possibilities for future research. Namely, we did not consider at all the questions of theoretical justification of viscopotential formulation and corresponding long wave models in the spirit of works by J. Bona, J.-C. Saut, D. Lannes, T. Colin [START_REF] Bona | Long wave approximations for water waves[END_REF]BCS02]. Other directions for future work consist in implementing new visco-potential framework in operational ocean and nearshore models.

Figure 1 .

 1 Figure 1. Conventional partition of the flow region into interior region and free surface, bottom boundary layers.

)

  This equation will be used in Section 5 to study the damping of linear progressive waves. Integral damping terms are reasonably well known in the context of KdV type equations. Various nonlocal corrections can be found, for example, in the following references[START_REF] Byatt-Smith | The effect of laminar viscosity on the solution of the undular bore[END_REF] Che68,[START_REF] Miles | Korteweg-de Vries equation modified by viscosity[END_REF][START_REF] Grimshaw | Damping of largeamplitude solitary waves[END_REF].A similar nonlocal KdV equation was already derived in [KM75]. They used a different scaling in boundary layer which resulted in dissipative term nonlocal in space. Later, Matsuuchi [Mat76] performed a comparison of numerical computations with their model equation against laboratory data. They showed that their model does not reproduce well the phase shift:

Figure 2 .

 2 Figure 2. Real and imaginary part of dispersion curve at t = 0. At the beginning, the nonlocal term has no effect. Thus, the real parts of the classical water wave problem and new set of equations are exactly superimposed on this figure. The imaginary part represents only local dissipation at this stage.

Figure 3 .

 3 Figure 3. Phase velocity at t = 1. Boundary layer effects start to be visible: the real part of the velocity slightly drops down and the straight lines of the imaginary part are deformed by the nonlocal term. Within graphical accuracy, the classical and their nonlocal counterparts are superimposed.

Figure 4 .

 4 Figure 4. Phase velocity at t = 2. Nonlocal term slows down long waves since the real part of the phase velocity decreases.

Figure 5 .

 5 Figure5. Phase velocity at t = 4. On Figure6we plot a zoom on long and moderate waves.

  2

Figure 7 .

 7 Figure 7. Limiting steady state of the phase velocity at t → +∞.

Figure 9 .

 9 Figure 9. Comparison among two dissipative and nondissipative Boussinesq equations. Snapshots of the free surface at t = 0.5

Figure 10 .Figure 11 .

 1011 Figure 10. Comparison among two dissipative and nondissipative Boussinesq equations. Snapshots of the free surface at t = 1.0

Figure 13 .

 13 Figure 13. Amplitude of the solitary wave as function of dimensionless time.

  Values of the parameters used in the numerical computations of the linear progressive waves amplitude. These values correspond to a typical Indian Ocean tsunami.Figure8. Amplitude of linear progressive waves as a function of time. Values of all parameters are given in Table1.

		parameter		definition		value
			ν		eddy viscosity		10 -3 m 2 s
			g		gravity acceleration	9.8 m s 2
			h			water depth		3600 m
			ℓ			wavelength		50 km
			k		wavenumber	= 2π ℓ m -1
	Table 1
					Amplitude of linear progressive waves
		1.02
		1
		0.98
	|A(t)|/|A(0)|	0.96
		0.94
		0.92
		0.9	0	1	2	3	4	5	6
						t, hours

  2νk 2 t = α 0 (1 -2νk 2 t) + O ν 2 k 4 t 2 .Figure 12. Comparison among two dissipative and nondissipative Boussinesq equations. Zoom on the soliton crest at t = 2.0

						t = 2.000
		2							No dissipation Local
		1.8							Nonlocal
		1.6
		1.4
		1.2
	η
		1
		0.8
		0.6
		0.4
		0.2
				0	0.2	0.4	0.6	0.8	1	1.2
							x
					Maximum wave amplitude as function of time
		2.01
		2
									Nonviscous model
		1.99
	max η(x,t)	1.97 1.98							Local terms
		1.96
		1.95
					Nonlocal equations
		1.94
		1.93	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8
							t

ℓ ∼ O(µ) .Inside the bottom boundary layer the normal gradient of the solenoidal part

Since the present problem is linear, we have effectively only two degrees of freedom but it is not important for our purposes.

To be precise, we plot real and imaginary parts of the dimensionless phase velocity which is defined as c p (k) := 1 √ gh ω(k) k .

In this section expression "classical" refers to complete water wave problem or its dispersion relation correspondingly.

Acknowledgment

The author of this paper would like to express his gratitude to professors Frédéric Dias and Jean-Claude Saut for very helpful discussions.