Antoine Varet

ENAC Nicolas Larrieu
email: nicolas.larrieu@enac.fr

Realistic network traffic profile generation: theory and practice

Keywords: Traffic generation, Network experimentation, Network performance analysis, Network measurement

Network engineers and designers need additional tools to generate network traffic in order to test and evaluate application performances or network provisioning for instance. In such a context, traffic characteristics are the very important part of the work. Indeed, it is quite easy to generate traffic but it is more difficult to produce traffic which can exhibit real characteristics such as the ones you can observe in the Internet. With the lack of adequate tools to generate data flows with "realistic behaviors" at the network or transport level, we needed to develop our tool entitled "SourcesOnOff". The emphasis of this article is on presenting this tool, how we implemented it and which methodology it follows to produce traffic with realistic characteristics. To do so, we chose to consider different stochastic processes in order to model the complexity of the different original traffics we wanted to replay. In our approach, we are able to consider several statistical laws and to combine their effects to model accurately the original behavior we analyzed in the real data. We then select the right parameters to consider as inputs for our SourcesOnOff tool. This approach gives really good traffic characteristics and, consequently, the generated traffic is really closed to reality as results presented at this end of this paper demonstrate it.

Introduction

Network engineers and designers need additional tools to generate network traffic in order to test and evaluate application performances or network provisioning for instance. In such a context, traffic characteristics are the very important part of the work. Indeed, it is quite easy to generate traffic but it is more difficult to produce traffic which can exhibit real characteristics such as the ones you can observe in the Internet.

With the lack of adequate tools to generate data flows with "realistic behaviors" at the network or transport level, we needed to develop our tool entitled "SourcesOnOff ". In this paper, "realistic behavior" means to generate traffic similar to the traffic a network administrator can capture on his network backbone, end-systems, etc... Previous research conducts us to design new methodologies and tools to generate network data traffic as close as possible to what we can find on Local Area Networks (LAN) and on the Internet. It can be necessary for instance to evaluate performances of new network entities (such as an embedded router [START_REF] Varet | Design and development of an embedded aeronautical router with security capabilities[END_REF]). In this specific context, we needed to face this router to generated traffic with characteristics as close as possible as the Internet traffic. However, depending on what you mean by "Internet network" and where you perform the study, measurements may completely differ: Internet cannot be solely characterized with a small set of parameters such as some mathematical distributions and additional factors. Indeed, there is not currently one unique mathematical modeling able to embrace the different characteristics and the complexity of the Internet traffic [START_REF] Olivier | Flow Level IP traffic characterization, France Télécom[END_REF]. Anyway, Internet traffics hopefully show some common properties and trends (detailed in the next paragraph) that helps to its modeling.

For instance, an Internet Service Provider (ISP) providing access for software engineering companies manage a different profile of data communications than an ISP for private individuals [START_REF] Gebert | Internet Access Traffic Measurement and Analysis[END_REF]. However, there are some common characteristics between both profiles. The goal of the tool we have developed is to handle most of Internet traffic profiles and to generate traffic flows by following these different Internet traffic properties.

The rest of this paper is structured as the following. The first section presents the mathematical model this tool is based on for traffic generation: the ON/OFF sources. We chose to consider different stochastic processes in order to model the complexity of the original traffic we want to replay. General approaches consider only one law for the On process and another law for the Off process. In our approach, we are able to consider several laws and to combine their effects to model accurately the original behavior we analyzed in the real data. We then select the right parameters to consider as inputs for our SourcesOnOff tool. This approach gives really good traffic characteristics and, consequently, the generated traffic is really closed to reality as results presented at this end of this paper demonstrate it. In the second section, we present why the other tools currently available were not usable in our case. The third section describes how our tool works plus some specific points the user should know about its implementation. The fourth section is dedicated to the validation of our tool with different methods in order to conclude if the generated traffic has the same properties than the original one. We have, in particular, investigated the traffic characteristics generated by our software and we demonstrate that it is able to generate network traffic with realistic properties such as these we can capture on the Internet. Different parameters have been considered to match generated traffic characteristics with the original traffic. First, we considered classical traffic parameters such as throughput, delay or losses but we have also investigated more advanced statistical parameters such as the correlation level of generated packets or their long range dependence (by computing the Hurst factor).

How to characterize an "Internet-like" profile?

Different studies were published with results related to Internet-like traffic characterization: the word "Internet" can cover very different profiles, but some common trends can be highlighted. Thus, we will examine two of them: high variability and self-similarity, respectively called the Noah and Joseph Effects.

High variability is characterized by an infinite mathematical variance and means that sudden discontinuous changes can always occur. Some mathematical distributions like Pareto and Weibull are heavy-tailed (i.e. the tail of the distribution is not exponentially bounded) and thus can be used to generate sets of values with high variances and also high-variability.

Different studies have shown that classical distributions such as the Gaussian and the Poisson distribution have failed to reproduce a LAN-like throughput [START_REF] Taqqu | On the Self-Similar Nature of Ethernet Traffic[END_REF] [START_REF] Olivier | Flow Level IP traffic characterization, France Télécom[END_REF]. One of them [START_REF] Olivier | Flow Level IP traffic characterization, France Télécom[END_REF] proposes to use Pareto's and Weibull's heavy-tailed distributions and proved that the generated traffic shows properties closer to real network traffic characteristics than other mathematical distributions.

Self-similarity is defined by a long-range dependence characteristic, which means there are bursts of traffic any time over a wide range of time scales. In other words, a small sub-range of values is "similar" to the whole range of values. W. Willinger found in [START_REF] Willinger | Self-Similarity Through High-Variability: Statistical Analysis of Ethernet LAN Traffic at the Source Level[END_REF] a relation between self-similarity and high variability for Ethernet Local Area Network (LAN) throughputs: in particular, he showed that the Noah and the Joseph effects are linked with a degree-1 polynomial relation. In other words, using ON/OFF sources with heavy-tailed distributions causes the traffic streams to be highly variable and, consequently, the aggregation of these streams to be also self-similar and highly variable.

ON/OFF sources to generate realistic throughputs

W. Willinger proposed in [START_REF] Willinger | Self-Similarity Through High-Variability: Statistical Analysis of Ethernet LAN Traffic at the Source Level[END_REF] to generate data following packet-train models, i.e. with strictly alternating, independent and identically distributed ON-and OFF-periods and proved that these ON/OFF sources generate data similar to experimental measurements on real networks.

A source corresponds to one network flow. The flow is associated with a source and destination couple and also with data, transmitted as a set of packets called a train of packets. In our case, we use the transport protocols TCP and UDP. The data stream is thus exchanged as a set of IP packets. The source is associated to a departure time and a duration time.

The departure time of any source is computed with the departure time of the preceding source plus a random duration. This randomness follows the user-defined "Doff distribution" (also named "distribution of inter-train durations"). The first source starts at the beginning of the process. Duration times should be similar for all sources. They should follow a user-defined random distribution called "Don distribution" (also named "distribution of train duration" or "distribution of flow duration"). This Don distribution is independent of the Doff distribution.

For feasibility constraints and to reduce the complexity of implementation, we use the Don distribution to generate random values, not for time duration in seconds but for quantity of transmission in bytes. Because of TCP congestion control mechanisms, source durations are correlated with quantities of transmitted data.

We can summarize the ON-OFF source generation algorithm we use in our tool by: They are issued from random number generators, which follow respectively the "Don" and the "Doff" distributions. These distributions should be heavy-tailed to ensure the self-similarity of the generated data throughput [START_REF] Willinger | Self-Similarity Through High-Variability: Statistical Analysis of Ethernet LAN Traffic at the Source Level[END_REF]. Distribution parameters are defined by the user on the command-line. Most of them provide "Internet-like" flow generators, some of them based on ON/OFF source generation processes. However, the aim of our research work was to face new network systems to real traffic. It means we wanted to perform our system evaluation in real time and do not want to evaluate a model of our system. Simulators cannot do this kind of real time experiments, only emulators can do it. However, it is complex to deal with most network emulators as they are complex to install, to configure with the real network environment, and consume as well many CPU resources. This is why we searched for and studied tools not to simulate or emulate but to generate a realistic traffic load on real networks.

Traffic replay tools

Some software enables the user to replay previously captured network traffics. This is done in two steps. Firstly, the user captures data on the networks he wants to reproduce. This can be done with tools called "sniffers" such as the reference tool Wireshark (Wireshark, 2013) and its background application "dumpcap" [START_REF] Dumpcap | Dump network traffic[END_REF], a command-line tool called up internally by Wireshark. Secondly, he needs to replay the captured traces, i.e. by retransmitting the sniffed packets in the same order, and separated with the same delays as these measured during the capture (cf. the "tcpreplay" (note 1) tool for instance). Harpoon (note 2) is another existing Open Source flow-level traffic generator, but this tool requires the user to define flow per flow the data weight he wants to transmit on the network. The automatic generation of the network profile is limited to constant and uniform distributions. Moreover, Harpoon seems to be not maintained since 2005. This is why we did not consider this tool for our experiments.

This traffic replay process has different advantages. The two steps may be done independently by two different users and at any time (as long as the first step is started before the second one). The capture may be filtered before being replayed. Replay parameters may be set up: for example traffic replay may be faster or slower than the capture speed, it may be altered with random losses or random transformations… However, replaying has also important drawbacks. Indeed, before replaying any network trace, the user must make or acquire the capture he wants to replay. Most of the time, this is a difficult task: user privacy issues restrain the administrators in allowing captures on the router they administer; often ISPs do not want to extract advanced statistics from their networks… The public data we have found on Internet are often so anonymized that they not longer contain the useful information we need to replay them. Moreover, if you continuously replay the same trace you will reproduce periodically the same "events" (throughput bursts, specific packet sequences or behaviors…). This periodicity might falsify some statistical analysis and deceive network systems. For instance, this type of traffic replay cannot be efficient if you want to assess performances of network mechanisms such as transport protocols given its objective is to act on traffic profile and by definition this traffic replay do not give any freedom to transport protocols to optimize traffic exchanges (such as congestion control mechanisms for instance).

Network throughput estimation tools

This is why different techniques and tools have been elaborated to load a network and evaluate its capacity. One of the existing tools we have studied is "iperf" [START_REF] Iperf | a modern alternative for measuring maximum TCP and UDP bandwidth performance[END_REF]. It can generate TCP or UDP flows to load the network. In TCP mode, an iperf client transmits to the server an infinite quantity of data through one TCP flow. After a user-defined duration, the iperf tool aborts the TCP connection and prints on the screen different statistics and the total quantity of data it succeeded in transmitting correctly.

A network capture shows often that mostly all the network resources are used by the TCP flows [START_REF] Gebert | Internet Access Traffic Measurement and Analysis[END_REF]. On one hand, the TCP protocol efficiently exploits the network and transmits data in an optimal time span. However, in LAN captures, external and unexpected events interact with the flow, they generate segment delays and losses impacting the TCP connection which becomes longer than the optimal time span expected. This drawback affects most of methodologies based on studies using only one TCP flow. We can cite, for example, when the tester uses an SSH or an FTP connection to transmit a big file and consequently to artificially load the network.

Figure 2c illustrates such a measure of throughput with iperf generating a TCP flow. We can see, from a first qualitative point of view, that the throughput generated by iperf is nearly constant and predictable, contrary to the throughput of an Internet capture (in figure 2a) and the throughput generated by our tool Sourceonoff (in figure 2b).

Iperf may be configured to transfer multiple TCP streams; in this case the program starts the different flows simultaneously at the beginning of the program. This behavior is not realistic: resources are shared fairly from the beginning, whereas in Internet later flows are penalized at their start by the flows already established.

Moreover, this tool "iperf" performs measurements on "long" TCP connections (about 10 second duration by default); this allows the TCP congestion control mechanism to adjust the data throughput efficiently. However, in most of ISP studies the TCP connections are very numerous and the majority is short (for instance, 81% of TCP connections carry less than 310 bytes, in study [START_REF] Iperf | a modern alternative for measuring maximum TCP and UDP bandwidth performance[END_REF]). Consequently, congestion control mechanisms do not have the time to optimize the connection.

In UDP mode, the client transmits data periodically: the period is computed depending on the user-defined quantity of data per second to transmit. However, for the same reason than previously, a single flow with constant throughput is not realistic [START_REF] Gebert | Internet Access Traffic Measurement and Analysis[END_REF] [START_REF] Olivier | Flow Level IP traffic characterization, France Télécom[END_REF].

Consequently, we have searched for other tools to generate data throughput: BWPing (BWPing, 2013), Ttcp [START_REF] Ttcp | test TCP/UDP performance[END_REF], NetPerf (NetPerf, 2013), NetPerfMeter [START_REF] Dreibholz | NetPerfMeter, A TCP/UDP/SCTP/DCCP Network Performance Meter Tool[END_REF], Ostinato (note 3)… All differ more or less on the set of supported protocols, on the maximum admissible throughputs and on the statistics printed on the screen for the user. However, all have the same drawbacks as explained above.

This is why, in the rest of this article, we chose the tool "iperf" as a comparison reference for our SourcesOnOff tool, because most measurement tools are based on the same assumptions and work in a similar manner: one or multiple identical and simultaneous TCP and UDP flows in order to load the network at its maximum available capacity.

Figure 2: Throughputs generated by iperf and our tool

The SourcesOnOff tool

We did not find in the state of the art we performed any tool to generate data flows based on exact ON/OFF sources. This is why we developed this one in C language and validated it with the Debian Operating System.

The tool we propose is free and Open Source, under the General Public License v3 (GPLv3 (GPL, 2013)). Source code can be downloaded at http://www.recherche.enac.fr/~avaret/sourcesonoff.

Initially developed for Debian-based systems, the 10 source files with 3,000 lines of code are compiled with the GNU C Compiler (gcc, with GCC (GCC, 2013)) and with the pedantic option. It is not POSIX-compatible because we use intensively the (long long) 64 bit-integer type and the (long double) float type for mathematical computations. Moreover, the tool uses certain Linux-specific functions (e.g. the epoll() system call) and so it should not be easily portable on non-Linux based systems, such as Microsoft Windows® or Mac OSX®. This software is linked with the real-time (-lrt) and the mathematical (-lm) libraries. The binary is between 84kB (32 bit version) and 111kB (64 bit version). Multiplying factors enables users to convert the randomly generated values into Bytes (for Don distributions), nanoseconds (for Doff distributions) and their multiples (kB, MB, GB, us, ms, s…). Distributions may be bounded by minimum and maximum user-defined values. For instance, the program may be limited to 1 GB per flow in order to avoid any traffic generation problems. It can be useful to set up heavy-tailed distributions. Indeed, most of them produce very high values such as 1e+44, but it is often problematic in practice: we cannot transmit a flow of 1e+44 bytes in our networks.

Our tool enforces user-defined minimum and maximum values for the generated random numbers, by increasing values lower than the minimum boundary and decreasing values greater than the maximum boundary.

An additional pseudo-distribution is available: the Constant distribution. In this case, all generated values are equal to a user-defined constant. This method enables the user to generate a more predicting behavior. For example, if the user calls the program for UDP flows with a Don distribution constant to 1kB and with a Doff distribution constant to 1 ms, the program will send each millisecond a datagram with 1024 bytes of data (note that this behavior is similar to iperf behavior in UDP mode).

A parameter called "turn number" enables the user to limit the number of flows: the program will start "turn" flows then wait for termination of all flows and quit. The default turn number is one million flows, in order to continuously create new flows for an time long enough (depending on the Doff values) to conclude most of our experiments. The special values of 1 flow for the turn parameter and a Don distribution constant to 1 GB is similar to transferring 1 gigabyte of random data through one TCP connection. In most cases, this overloads the network capacity of the path until all the data are transmitted.

Generation of the network flows

First of all, different sets of Don and Doff random values are generated. They are then used for data communications. Command-line parameters specify, for each set of (Don; Doff) values, the different socket properties: are they transmitters and so who is the destination host? Or is it a bottomless pit? Do we use UDP or TCP protocol? What is the port number to use? Does IPv4 or IPv6 be used with? Does a set of sources be delayed before starting (to let another set of sources take some advance)? Does the set of sources be stopped after a delay, even if there are still active sources? How many maximum sources should be created? Does the random generator be initialized with a specific seed (i.e. to reproduce two exact source generations)?

All these parameters enable the user to refine the tool behavior. Different sets of sources may run simultaneously, each set is associated with independent Don and Doff distributions and parameters. For instance, a set of sources may be defined to load the network with TCP streams and another set to generate UDP datagrams. The process forks itself to provide TCP and UDP data simultaneously. Each task uses the epoll() primitive to perform the management of thousand of flows if necessary. This primitive reduces the portability of the software to only Linux systems and should be replaced with the POSIX-compatible select() primitive in a future tool release.

Contrary to programs like iperf, this tool is only intended to generate data. It was not developed to provide itself advanced statistics on the network: we use our tool to generate background traffic and additional passive and active measurement tools to evaluate network performances and collect statistics (note 4).

Statistic profile extraction from a real traffic trace

Before experimenting traffic generation, we have to define what kind of traffic we want to generate. This step may be done arbitrarily by choosing appropriate values or by selecting values from existing studies like [START_REF] Gebert | Internet Access Traffic Measurement and Analysis[END_REF] [START_REF] Olivier | Flow Level IP traffic characterization, France Télécom[END_REF].

We chose an alternative solution: given that we belong to a research entity, we have easily access to real data. Thus, we asked to our local network administrator to capture the entire incoming and outgoing traffic, generated by people from the university (students, administrative people, teachers and researchers) on our local area network. The data were captured on the firewall protecting the access link between our LAN and other networks: REMIP (REMIP, 2013) and RENATER [START_REF] Renater | [END_REF]. REMIP and RENATER networks are our links to the Internet. Thus, we captured all data from our LAN to Internet and all associated responses.

We captured a set of 9 millions of IPv4 packets, mostly TCP data (97.7% of TCP, 2.2% of UDP and 0.1% of ICMP), during 10 hours between 8:00AM and 6:00PM, Tuesday the 29th of January, 2013. We expected to generate simultaneously TCP sources and UDP sources with our tool SourcesOnOff. However, UDP datagrams are negligible in our case (less than 2.5% of our traffic in packet numbers, less than 0.5% of our traffic in bytes), so we replayed only TCP sources. We will finally just model and reproduce the TCP traffic in our study case (section 4 will explain it). We used the tool tcpdump [START_REF] Tcpdump | Tcpdump & LibPcap public repository[END_REF] to capture the raw packets and then the tool ipsumdump (IPsumDump, 2013) to ensure the anonymity. We finally conducted a complete statistical analysis with bash and R scripts. These tools enabled us to retrieve different statistical characteristics of the captured data. This process is detailed in the next subsection.

Statistic profile extraction process

Previous research work (such as [START_REF] Gebert | Internet Access Traffic Measurement and Analysis[END_REF] [START_REF] Olivier | Flow Level IP traffic characterization, France Télécom[END_REF] (Leland et al., 1997) [START_REF] Taqqu | On the Self-Similar Nature of Ethernet Traffic[END_REF] for instance) have demonstrated that one unique statistical law cannot figure out all the complexity of an original Internet traffic trace. This is why we chose to decompose and to model the Internet traffic trace we want to replay by using several different statistical laws. Thus, we need, firstly, a decomposition algorithm and, secondly, a distance criterion to evaluate the differences between real original data and data generated by our tool. This distance criterion will help us to select the best statistic laws to decompose original data. These two parts of the process are detailed in the two next subsections.

Traffic trace decomposition

We developed an algorithm able to detect a lack of continuity in any data we want to characterize. This algorithm is mainly based on the quantmod tool developed by Jeffrey A. Ryan (note 5). We explain the detection results on the example provided by figure 4. The top part of the figure represents the original data (before applying our decomposition algorithm) and the bottom part of the figure represents the decomposed data. On the left side, we represent the ECDF of the data and on the right side we compare the original data to a specific statistic profile (we chose to show only the Weibull function comparison in this figure) thanks to QQPlot diagrams (Becker et al., 1998). We can note that before decomposition the statistic profile does not fit the whole set of data (cf. the top QQPlot diagram where a linear tendency cannot be seen). Indeed, for this trace profile, a Dirac distribution creates, at its origin, a discontinuity in the original data. By using our detection algorithm we are able to generate the same original statistical process by considering two statistical sub processes. The first one is the Dirac function that creates this discontinuity and the second one is the relevant function that fits the data without the Dirac function's events. Thus, in the bottom part of figure 4, we consider finally a Weibull distribution for the remaining data and the reader can note that it perfectly fits its statistic profile (cf. the bottom QQPlot diagram where a linear tendency can be seen).

To characterize this resulting process of the remaining data, we generated distributions of different types (Weibull, Pareto, Exponential, Gaussian…) and with different parameters. For instance, the script we wrote to optimize Weibull shape parameter tests all Weibull values between 0.01 and 10.00 with a step of 0.01. For each of the 10.000 candidate values, we measured the distance between the network trace values and an equivalent number of values following the Weibull distribution with the candidate shape. We then select the shape value with the minimum distance, in the 10.000 computed correlations.

It is now mandatory to explain how this distance assessment is computed. For this purpose, we used the Bayesian Information Criterion which is introduced in the next subsection. We are thus capable of combining the different values of the traffic to generate by considering in the SourcesOnOff tool those different distributions. The next step of the generation process is to quantify the statistical distance that exists between our original data and the data generated by our tool. To compute this distance we introduce the Bayesian Information Criterion [START_REF] Schwarz | Estimating the dimension of a model[END_REF]. This criterion is computed as BIC = k * ln(n) -2 * ln(L), where:

-n is the size of analyzed data; -L is the likelihood of the model (Weibull, Pareto, Exponential…) regarding the different original data; -k is the total number of estimated parameters.

The final goal of this comparison is to select the smallest BIC (minimum BIC value is -∞) according to the different candidate statistic profile (Weibull, Pareto, Exponential, Gaussian…). Thus, the tool can conclude that the selected model (which might be a composition of different distributions) is the closest from the original data.

Validation of the SourceOnOff tool

The objective of this section is to validate the traffic profile generated by our tool. We could have presented complex and realistic topologies, but the different experiments we performed have not given interesting results and network behaviors. Moreover, the goal of this section is to analyze deeply the traffic generation process provided by our tool. By introducing additional complex topologies, the statistical analysis of the generated traffic would have been more complex given that we would have had difficulties to link a parameter variation with a specific cause: our tool or the complex topology. This is why we chose a simplified experimental topology with only two hosts and one router, as presented in figure 5. The SourcesOnOff tool has been deployed on the different hosts: the sender part on the transmitter host and the receiver part on the receiver host. This subsection details the deployment of an expected traffic profile on the different end-systems. For this experiment, we use two Linux Debian network hosts called H1 and H2, running on a mono-core 1.2 GHz with 1 GB of RAM. They are connected to a router RT running a quad-core Intel processor @3.8 GHz with 4 GB of RAM. The links between hosts are Ethernet RJ45 connections manually configured to 10baseTx-HD, in order to easily limit the link capacities without any software solution and thus to capture all the traffic without any loss. The experimental network topology has been already presented in figure 5. H1 runs the sender part of our tool and H2 runs the receiver part.

After the beginning of the generation, the traffic profile tends quickly to a permanent and stable mode. At this moment, the processor consumption reported by system monitoring tools such as top & htop is around 5% CPU time on the transmitter H1 and 10% on the receiver H2. The virtual memory allocated by the Linux Operating System for the tool is between 2 and 16MB on our hosts, depending on the number of sources we wants to start (16 MB when infinite number of sources). These measures remain stable until the program automatically finished or was stopped by the user.

We captured different network traces: between 10 minutes and 10 hours. We will describe in the next subsections only one of them, because other captures concluded to the same results. The following subsections will analyze the captured data on the router during one hour, searching for the properties which characterize an Internet-like traffic profile.

Verification of generated traffic correctness

The validation of our tool needs to analyze its generated traffic. This is why this section analyzes the original captured traffic and explains how generated traffic is really close to the original network traffic. Thus, we performed different quantitative and qualitative verifications on the generated traffic to answer the following questions: are the requested parameters for the Don and Doff distributions correctly applied by the tool, are the random values conforms to the "original" real capture, and finally, is the generated traffic conform to the original one?

Statistic profile detection

By applying the algorithm described in section 3.3, the tool is able to detect several different statistic profiles in the same original data as described in the Figure 6. The left side describes the Don process. We can graphically see that both original and generated distributions are well fitted by the Weibull function (the "generated" black curve fits perfectly the "original" red one). For the Doff distribution it is a more complex process that we have modeled based on the function composition algorithm. Indeed, it chose a Weibull function plus a Dirac function (as plotted by the black curve). In the following sub section, we are going to validate how this function composition fits very well the original statistical process.

Quantile-quantile plots

As qualitative estimators, we used quantile-quantile plots diagrams (Becker et al., 1978) (also called "QQPlots") to represent on the same diagram the measured values on the real network (i.e. the original data) and the measured values on the experimental network (i.e. the SourcesOnOff generated data).

QQPlot diagrams sort independently X and Y values, and then represent points to the sorted coordinates. When the X and Y series are correlated, a linear tendency is visible on the diagram. On our QQPlot diagrams, we will represent a diagonal blue line of equation y=x, in order to represent a perfect correlation between x and y values.

Firstly, we studied the Doff values, i.e. the inter-train durations. The set of values used for the X data is the set of durations between two consecutive TCP connections observed in the capture in the real network. The set of values used for the Y data is the set of durations between two consecutive TCP connections measured in the traffic generated by our tool SourcesOnOff. In other words, the Y values are the durations our tool waited between starting two TCP connections.

We can see on figure 7 that the values are very well correlated: most of the 9,700 points drawn on this figure are near the diagonal blue line of equation y=x. This means most values measured in the real network were successfully reproduced in our experimental network. The correlation factor is 99.8% in our case for Doff values. A few points with high values are not well correlated: they come from the tail of the Weibull distribution densities.

Secondly, we performed the same analysis with the Don values, i.e. the measurement of TCP connection durations (exchanged bytes between the different hosts). Figure 8 shows the QQPlot diagram. Observed values are well correlated with expected values. The correlation factor is 97.9% in our case for Don values.

We can conclude that our tool shows very good trends for the generated traffic compared to the original one. However, it is necessary to analyze more accurately how the generated traffic is close to original. This is why autocorrelation function is described in the next subsection.

Autocorrelation checking

We computed the autocorrelation on the throughput generated by our tool SourcesOnOff. We compared this autocorrelation function with the autocorrelation on the Internet throughput we captured on our real network. Each time, we computed the autocorrelation with a maximum lag of 15,000 values.

We can see on figure 9 that autocorrelation functions show both heavy tails: when the degree of correlation increases, the value of autocorrelation stays around 0 but does not converge to 0. Instead, the autocorrelation oscillates around 0 with bursts upper the limit of 5% materialized by the blue dotted lines. This is a qualitative proof of the long-term memory feature on the measures. However, we use additional tools (BIC distance and Hurst factor) to check quantitatively these results in the following section.

BIC distance

Table 1 shows BIC distance analysis. We can note that original and generated traffics are the closest when the statistic model combines both Weibull and Dirac profiles. The qualitative similarity between original and generated traffics provided by QQPlot is thus confirmed by BIC distance computations for both Don and Doff profiles. This result validates that a combination of both Weibull and Dirac distributions is the right statistical profile to consider with our original data.

Hurst exponent computation

We wanted also to compute quantitatively the long-term memory of the generated traffic compared to the original one. This can be showed with the Hurst exponent computation: this value is an indicator of long-term memory. In our case, we used the Wavelet-based joint estimator of the Hurst exponent, as described by D. Veitch and P. Abry in [START_REF] Veitch | A Wavelet Based Joint Estimator of the Parameters of Long-Range Dependence[END_REF] and implemented in C in [START_REF] Roughan | Program in C to estimate on-line the hurst component[END_REF]. This estimation algorithm is the one used in the tool LDEstimate (LDEstimate, 2012) to study the Long-Range Dependence. In table 2, we compare the Hurst exponent estimation for the throughput generated by our tool with the throughput of the real capture at different time scales (from 100 us to 10 s).

Conforming to [START_REF] Veitch | A Wavelet Based Joint Estimator of the Parameters of Long-Range Dependence[END_REF], data (see table 2) exhibit Long Range Dependence (LRD). Indeed, data show LRD when 0.5 < H < 1. Moreover, the dependence is stronger when H is closer to 1. Thus, we can conclude than in our case, we have a strong Long Range Dependence in both series of throughput, whatever the time window is defined. We can also conclude that both generated and original traffics exhibit the same LRD level (considering a 12 % error interval).

SourcesOnOff usages

In this section, we address specific usages that we can imagine for this tool. For instance, it may be useful to apply transformations on the traffic profile we want to generate. Thus, let us imagine an example where we have a real data and its extracted statistic profile to be equivalent to a load of 100 users on the network. Assuming the hypothesis "the throughput is proportional with the number of users" and "the inter-flow durations are inversely proportional to the number of users", we can do network provisioning with some simple mathematical operations on the raw parameters of the profile. It is, for us, extremely harder to do the same operations with a capture of raw packets that we would like to replay.

In this section, we will consider an example where we want "to emulate the network intensity in the case we have now 150 users (rather the 100 users observed in the real data) but the same systems". The goal is to be able to predict the final throughput Y (Mbps) that we need to generate with our tool to represent the 150 users.

If the parameters of the Don and the Doff distributions are not adequate, the network throughput may differ greatly from expected. We can do an analogy with a sinusoidal signal. As illustrated in figure 10, if the signal is too low, the capacity will be underused, while if the signal is too much amplified, it will saturate the channel and measures will differ from a sinus.

Figure 10: Explanation of tuning with the sinus function With a high-capacity network infrastructure, high values for inter-train durations (Doff) and low values for quantity of data per flow (Don), the network throughput will mostly be empty (~0% load), with bursts for each source starting. On the contrary, with a low-capacity network infrastructure, high values for Don and low values for Doff, the network will be soon saturated (~100% load). Then, the transmitter system may fail to create new sources, due to not enough network resources to establish new connections, until old connections finishes.

Figure 11 illustrates the three cases of tuning: the blue line at left represents a low throughput generated by SourcesOnOff, the black line at the center of the figure a medium throughput and the blue line at right a high overhead where the network path is mostly saturated. Note that each of the three cases provides some information: high throughputs enable the tester to estimate the network path capacity, low throughputs to evaluate system behavior in front of burst and medium throughputs to evaluate system behavior to a regular and long stress.

We can expect an average throughput of X Mbps = Don.average * (1/ Doff.average). In order to realize an average throughput of Y Mbps instead of the current average of X Mbps, we have to adjust proportionally the Doff average value: Doff_new_mean = (X/Y)*Doff_old_mean. In our opinion, the Don values correspond to the behavior of the users and the Doff values is inversely proportional to the number of users. So changing the Don parameters would make the flows unrealistic whereas the Doff average value would reduce the number of simultaneous flows.

For some specific statistical profiles, the mean is not a characterization parameter of the distributions. For instance, this is the case for the Weibull distributions. These distributions use mostly a "scale" parameter which is proportional with the mean of the distribution. Thus, with Weibull distributions, we should apply the following equation: Doff_scale_new = (Y/X)*Doff_scale_old.

If we come back to the pratical example we are considering in this section, we should do the same operation to reach our objective "emulate the network with 150 users instead of the 100 users". This is achieved by dividing the Doff scale by a factor of 1.5 and thus we use a scale factor of 193 ms for the Doff Weibull distribution. This simple example gives a overview on how this SourcesOnOff tool can be easily tuned in order to manage different performance goals.

Conclusion and future work

In this paper, we have introduced a methodology to generate network traffic with realistic characteristics. A tool has been developed, based on the application of ON/OFF sources with different statistical profiles. Parameters of the distributions can be defined by the user or extracted from real traffic analysis. We have completed this paper with a validation of both the traffic generation methodology and the SourcesOnOff tool. Different experiments have been conducted. All of them validate that our tool is able to generate traffic with the same characteristics than real ones.

The tool is freely available and may be utilized for a wide variety of network traffic profiles. We hope users can appreciate the wide range of applications where SourceOnOff can be utilized. We can cite, for instance, network provisioning, application testing or network mechanism enhancement. For this last field of application, this tool can be useful to generate random data in order to observe bufferbloats, to stress an intermediary network system with a big number of simultaneous flows or to generate network traffics to observe TCP race conditions. However, additional work would have to be conducted by considering different traffic traces collected in different points of presence of the Internet. We could imagine to have an access in the future to traffic captured not only on the firewall of a LAN but also on a core router of the Internet. By putting our SourceOnOff tool freely available to the research community, we hope we can provide a first level of contribution to improve characterization results for Internet traffic understanding. Furthermore, this tool is currently unable to generate traffic different than TCP and UDP protocols. In the future, the tool may be completed with support for ICMP protocol or additional transport protocols such as SCTP (Stream Control Transmission Protocol) [START_REF] Stewart | Stream Control Transmission Protocol[END_REF] for instance. Moreover, the tool may support other statistical distribution profiles and may provide additional statistics for the users. The validation of this tool being completed, we can now take into account more complex network topologies (cloud computing applications for instance) and distribute different SourcesOnOff sender and receiver agents among them for the future experiments we plan to realize.

  For example, if we have the following sequence of sources <A, B, C, …> and if Doff values, randomly generated, begin with the sequence <A.Doff=3, B.Doff=3, C.Doff=5, …> and Don values with <A.Don=10, B.Don=3, C.Don=4, …>, then figure 1 illustrates which throughput we can expect on the network.Internet flows are mainly TCP and UDP flows. From an application point-of-view, the TCP protocol uses congestion avoidance mechanisms and then spreads the data transmission out in time: the stream is divided into segments and each segment is sent when the previous segment has been correctly received. Flow durations depend of many unpredictable factors in real systems where other flows share the same resources. Moreover, the UDP protocol divides the data into multiple datagrams and sends all datagrams as fast as the lower layers can, with possible losses. The proportion of TCP vs. UDP is a parameter of the Internet profile.

Figure 1 :

 1 Figure 1: Emission process of ON/OFF traffic sources 2. State of the art of traffic generators 2.1 Network simulators Different tools are available to simulate networks and their behaviors, with integrated ON/OFF sources. We can cite, for instance, NS-2 (NS-2, 2013), OpNet (OpNet, 2013) or OmNet++ (OmNet++, 2013).

Figure 3 :

 3 Figure 3: SourcesOnOff software design (grey box view) The source code is divided into 3 different tasks represented in figure 3: the random value generation task, the network flow generation task and a third task to link the two previous together (common declarations and structures, command-line parsing…).

3. 1

 1 Generation of the random values ON/OFF sources are based on the generation of random values following well-determined distributions. The following distributions are currently implemented in the program. The Uniform distribution is based on the drand48() and the random() Linux functions, with corrections to ensure an exact uniformity on any range. The Normal/Gaussian distribution is computed with the help of the Box-Muller transformation. We use the Knuth's algorithm for the Poisson distribution. The Pareto, Weibull and Exponential distributions use internally a transformation on the uniform distribution. The generation process of all these distributions have been statistically validated with the R statistical software (R, 2013) and its comparison function qqplot().

Figure 4 :

 4 Figure 4: Example of traffic decomposition (ECDF [left] and QQPlot [right] functions for real data [top] and Weibull distribution [bottom]) 3.3.3 BIC (Bayesian Information Criterion) distance assessment

Figure 5 :

 5 Figure 5: Experimental network topology 4.1 Tool deployment

Figure 6 :

 6 Figure 6: Fitting of original (red curve) and generated (black curve) traffic distributions (left side Don and right side Doff distributions) 4.2.2 Qualitative analysis

Figure 7 :

 7 Figure 7: QQPlot for Doff values (real traffic vs. generated traffic)

Figure 8 :

 8 Figure 8: QQPlot for Don values (real traffic vs. generated traffic) 4.2.3 Quantitative analysis

Figure 9 :

 9 Figure 9: Generated (top) vs. real (bottom) throughput autocorrelations

Figure 11 :

 11 Figure 11: Explanation of tuning impact on SourcesOnOff 5. Conclusion and future work

Table 1 :

 1 BIC distance between original and generated traffic

	Statistic model	Doff BIC distance	Don BIC distance
	Weibull (+ Dirac)	0.01558	0.02527
	Pareto (+ Dirac)	0.0209	0.03834
	Exponential (+ Dirac)	0.08892	0.06222
	Gaussian (+ Dirac)	0.1027	0.07535

Table 2 :

 2 Hurst exponent estimations

	Sample duration for	Hurst Exponent estimation	H exponent for the	Ratio between generated
	throughput	(H)	real throughput	traffic and real traffic
	measure-ments	for our tool
		SourcesOnOff
	100 us	0.88	0.83	7 %
	10 ms	0.97	0.95	2 %
	100 ms	1.00	1.00	0 %
	1 second	1.00	0.88	12 %
	10 s	1.00	1.00	0 %

Acknowledgements

We would like to thank Benoit Sainct (MS student) for his work in implementing some of the statistical tools used for traffic characterization.

Copyrights

Copyright for this article is retained by the author(s), with first publication rights granted to the journal.

This is an open-access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/3.0/).


C Cherubini
email: claudia.cherubini@gmail.com

C I Giasi

N Pastore

Bench scale laboratory tests to analyze non-linear flow in fractured media

The knowledge of flow phenomena in fractured rocks is very important for groundwater resources management in hydrogeological engineering.

A critical emerging issue for fractured aquifers is the validity of the Darcian-type "local cubic law", which assumes a linear relationship between flow rate and pressure gradient to accurately describe flow patterns.

Experimental data obtained under controlled conditions such as in a laboratory increase our understanding of the fundamental physics of fracture flow and allow us to investigate the presence of non-linear flow inside fractures that generates a substantial deviation from Darcy's law.

In this study the presence of non-linear flow in a fractured rock formation has been analyzed at bench scale in laboratory tests. The effects of non-linearity in flow have been investigated by analyzing hydraulic tests on an artificially created fractured rock sample of parallelepiped (0.60 × 0.40 × 0.8 m) shape.

The volumes of water passing through different paths across the fractured sample for various hydraulic head differences have been measured, and the results of the experiments have been reported as specific flow rate vs. head gradient. The experimental results closely match the Forchheimer equation and describe a strong inertial regime. The results of the test have been interpreted by means of numerical simulations. For each pair of ports, several steady-state simulations have been carried out varying the hydraulic head difference between the inlet and outlet ports. The estimated linear and non-linear Forchheimer coefficients have been correlated to each other and respectively to the tortuosity of the flow paths.

A correlation among the linear and non-linear Forchheimer coefficients is evident. Moreover, a tortuosity factor that influences flow dynamics has been determined. C. Cherubini et al.: Bench scale laboratory tests to analyze non-linear flow in fractured media   

Introduction

The purpose of this paper is to experimentally investigate the behavior of high velocity flow regimes in a fractured network at bench scale.

High velocity flow dynamics can have significant impact in diverse fields such as radioactive waste disposal, geothermal engineering, environmental risk assessment and remediation, reservoir engineering, and groundwater hydrology [START_REF] Cherubini | Fluid flow modeling of a coastal fractured karstic aquifer by means of a lumped parameter approach[END_REF].

In most studies examining hydrodynamic processes in saturated porous and fractured media, it is assumed that flow is described by Darcy's law, which expresses a linear relationship between pressure gradient and flow rate. Darcy's law has been demonstrated to be valid at low flow regimes (Re ≪ 1). For Re > 1 a non-linear flow behavior is likely to occur. Ing and Xiaoyan (2002) have shown how non-Darcian flow has a significant impact on consolidation rate in geotechnics. [START_REF] Basak | Effect of non-Darcy flow on computing seawater intrusion lenghts[END_REF] have demonstrated that seawater intrusion length increases when flow deviates from Darcian linearity.

Exposure and risk assessment of chemical pollution based on the applicability of Darcy's law are sometimes inappropriate when applied to fractured aquifers because of the Published by Copernicus Publications on behalf of the European Geosciences Union.

C. Cherubini et al.: Bench scale laboratory tests to analyze non-linear flow in fractured media

possible occurrence of non-laminar flow regimes. In fact, considerable evidence exists to refute their reliability in assessing solute-migration rates and for determining downgradient concentrations [START_REF] Field | Risk Assessment Methodology for Karst Aquifers (2) Solute transport Modeling[END_REF]. In solute transport, non-Darcian flow might give rise to tailing in contaminant breakthrough curves that show a non-Fickian behavior [START_REF] Boutt | Trapping zones: The effect of fracture roughness on the directional anisotropy of fluid flow and colloid transport in a single fracture[END_REF][START_REF] Cardenas | Navier-Stokes flow and transport simulations using real fractures shows heavy tailing due to eddies[END_REF].

A non-Fickian behavior is also valid in case of highconcentration brine transport, where the assumption is that the linear Darcy's law holds as shown in [START_REF] Watson | On the validity of darcy's law for stable high-concentration displacements in granular porous media[END_REF]. On the other hand, [START_REF] Tenchine | Density contrast effects on tracer dispersion in variable aperture fractures[END_REF] carried out density driven flow simulations in a rough-walled natural fracture extracted from a limestone quarry and observed measurable non-linearity between the velocity growth rate and the velocity, indicating that the Navier-Stokes formulation must be used.

The mathematical representation of fluid dynamics in fractured rock aquifers is of a great concern for environmental and petroleum engineering and in geological sciences.

The local cubic law adapts Darcy's law to flow through fractures under the assumption of ideal fractures with flat, smooth and parallel walls and infinite lengths, together with laminar flow, incompressible fluid and confined aquifer configuration.

Real rock fractures instead are characterized by rough walls and variable surfaces, geometry and apertures. The presence of asperities and obstructions or sharp changes in fracture profile is the reason for microscopic inertial phenomena that cause an extra macroscopic pressure loss which deviates flow from linearity.

Roughness has a large influence on fluid flow and transport through tight, rough-walled fractures [START_REF] Boutt | Trapping zones: The effect of fracture roughness on the directional anisotropy of fluid flow and colloid transport in a single fracture[END_REF] where non-Darcian flow is particularly easy to occur (e.g. [START_REF] Lomize | Filtratsia v treshchinovatykh porodakh, Seepage in Jointed Rocks, Gosudarstvennoe Energeticheskoe Izdatel'stvo[END_REF][START_REF] Louis | A study of groundwater flow in jointed rock and its influence on the stability of rock masses[END_REF][START_REF] Qian | Experimental study of turbulent unconfined groundwater flow in a single fracture[END_REF][START_REF] Qian | Experimental evidence of scale dependent hydraulic conductivity for fully developed turbulent flow in a single fracture[END_REF][START_REF] Qian | Experimental study of solute transport under non-Darcian flow in a single gracture[END_REF][START_REF] Wen | Non-Darcian flow in a single confined vertical fracture tward a well[END_REF].

Tortuosity, which usually characterizes the ratio of the effective path length connecting two locations in porous media to the geometric distance [START_REF] Tenchine | Density contrast effects on tracer dispersion in variable aperture fractures[END_REF], has been found to significantly affect fluid flow in fractured media under certain conditions [START_REF] Tsang | The effect of tortuosity on fluid flow through a single fracture[END_REF][START_REF] Wang | Hydrologic mechanisms governing fluid flow in apartially saturated fractured porous medium[END_REF]. [START_REF] Yeo | Applicable range of the Reynolds equation for fluid flow in a rock fracture[END_REF] identified a criterion parameter function of the roughness and tortuosity for the applicability of the Reynolds equation to fluid flow in rock fractures. [START_REF] Louis | Etude des écoulements d'eau dans les roches fissurées et de leurs influences sur la stabilité des massifs rocheux[END_REF] empirically defined five steady-state flow regimes, from smooth laminar to fully developed turbulent, depending on various degrees of relative roughness, fracture aperture, hydraulic head gradient in the fracture plane and kinematic viscosity.

Fully developed turbulent flow conditions are more likely to occur when velocities are high, i.e. in karstic conduits, whereas, in presence of lower velocities, non-linear laminar flow happens when inertial effects become important, e.g. due to roughness, restrictions, presence of contact points between the fracture walls and non-rectilinear profile.

Different from conduits, the flow regime prevailing between rough-walled surfaces does not immediately switch from Darcy flow to turbulent flow as velocity is increased or vice versa [START_REF] Kohl | Observation and simulation of non-darcian flow transients in fractured rock[END_REF], but transitional regimes can be distinguished.

On the basis of a recent classification by [START_REF] Skjetne | High-velocity flow in porous media, analytical, numerical and experimental studies[END_REF], three different principal flow regimes can be distinguished:

(1) linear laminar flow, (2) non-linear rough laminar flow (weak or strong inertia) and (3) turbulent flow.

A fracture can be envisioned as two rough surfaces in contact. Cross sectional solid areas representing asperities in contact are similar to the grains of porous media. It is therefore possible to apply the general equations describing flow and transport in porous as well as in fractured media. The void space of the fracture is interconnected in one surface, so the fracture is best described as a two-dimensional porous medium (National Research Council, 1996).

In the literature different laws are reported that account for the non-linear relation between velocity and pressure gradient.

The weak inertia equation is a cubic extension of Darcy's law and describes pressure loss versus flow rate for low flow rates:

-∇p = µ k • v + γρ 2 µ • v 3 (1)
where

p (M L -1 T -2 ) is the pressure, k (L 2 ) is the permeabil- ity, µ (M L -1 T -1 ) is the viscosity, ρ (M L -3 ) is the density, v (L T -1
) is the velocity, and γ (L) is called the weak inertia factor. This law was first shown numerically [START_REF] Barrére | Modélisation des écoulements de Stokes et Navier Stokes en milieu poreux[END_REF]) and then derived theoretically for homogeneous isotropic media [START_REF] Mei | The effect of weak inertia on flow through a porous medium[END_REF].

In case of higher Reynolds numbers (Re ≫ 1), the pressure losses pass from a weak inertial to a strong inertial regime, described by the Forchheimer equation [START_REF] Forchheimer | Wasserbewegung durch Boden[END_REF], given by:

-∇p = µ k • v + ρβ • v 2 (2)
where β is called the inertial resistance (coefficient), or non-Darcy coefficient.

Forchheimer proposes a further expression in which the pressure gradient can be described by a third order polynomial, whereas [START_REF] Hassanizadeh | High velocity flow in porous media[END_REF] propose an extension that includes the effects of a transient flow regime.

Another commonly used non-Darcian flow equation is the Izbash or power-law equation [START_REF] Izbash | O filtracii V Kropnozernstom Materiale[END_REF]:

-∇p = λ • v n (3)
where λ and n are two constant coefficients and 1 ≤ n ≤ 2. If n = 1 the flow is Darcian, whereas if n = 2 the flow is fully developed turbulent. The Forchheimer and Izbash equations are equivalent when the linear term in the Forchheimer equation is equal to zero and n = 2. The Forchheimer equation can be considered as an expansion of the linear Darcy's law because it includes both the viscous and inertial effects.

On the other hand, power law functions such as the Izbash equations are more appropriate to model post-linear non-Darcian flows that might be caused predominately by turbulent effects [START_REF] Wen | An analytical solution for non-Darcian flow in a confined aquifer using the power law function[END_REF].

Experimental evidence has shown that both the Forchheimer and the Izbash equations are equally possible, depending on field-specific conditions [START_REF] Wen | Non-Darcian flow in a single confined vertical fracture tward a well[END_REF].

A general Darcian-like relationship can be used [START_REF] Chin | Nonlinear Flow in Karst Formations[END_REF] to describe all the mentioned flow regimes and to account for non-linearities in the relationship between hydraulic head gradient and flow velocity:

v = -K eff (∇h) • ∇h (4)
where K eff (L T -1 ) is the effective hydraulic conductivity, and h = p/ρg (L) is the total hydraulic head. For instance, according to Forchheimer's law, the effective hydraulic conductivity can be written as [START_REF] Cakmak | Analysis of nonlinear darcy -forchheimer flows in porous media[END_REF]:

K eff = 2 a + a 2 + 4b ∇h (5)
where a (T L -1 ) and b (T 2 L -2 ) are the linear and inertial coefficients, respectively, in terms of hydraulic head and can be expressed according to the previous expression as:

a = µ ρgk ; b = β g . (6)
In the same way the effective fracture transmissivity for a discontinuity can be defined, and the relationship between specific flow rate q (L 2 T -1 ) and hydraulic head gradient can be written as:

⇀ q = 2 a f + a 2 f + 4b f ∇h ∇h (7)
where a f and b f are related to a and b:

a f = a w ; b f = b w 2 (8)
where w represents fracture aperture. Flow regimes and non-linear behavior of fluid flow through fractures have been investigated empirically [START_REF] Lomize | Filtratsia v treshchinovatykh porodakh, Seepage in Jointed Rocks, Gosudarstvennoe Energeticheskoe Izdatel'stvo[END_REF][START_REF] Louis | A study of groundwater flow in jointed rock and its influence on the stability of rock masses[END_REF], experimentally [START_REF] Witherspoon | Validity of cubic law for fluid flow in a deformable rock fracture[END_REF][START_REF] Qian | Experimental study of turbulent unconfined groundwater flow in a single fracture[END_REF][START_REF] Chin | Nonlinear Flow in Karst Formations[END_REF]numerically (Zimmerman et al., 2004;[START_REF] Kolditz | Non-linear flow in fractured rock[END_REF][START_REF] Brush | Fluid flow in synthetic roughwalled fracture: Navier-Stokes, Stokes, and local cubic law simulations[END_REF].

Zimmerman et al. ( 2004) studied non-linear flow regimes both with experimental and numerical simulations of Navier-Stokes equations and suggested the critical Reynolds number of 10 for practical purposes. 2.44 Hydraulic conductivity (m s -1 ) 1.63e-8 [START_REF] Witherspoon | Validity of cubic law for fluid flow in a deformable rock fracture[END_REF] conducted experiments in rock fractures to study the hydraulic behavior. They demonstrated that if fractures are rough or flow velocities are high, Darcy's law is not applicable. [START_REF] Chin | Nonlinear Flow in Karst Formations[END_REF] investigated the variation of effective hydraulic conductivity as a function of specific discharge in several 0.2 m and 0.3 m cubes of Key Largo Limestone. The experimental results closely match the Forchheimer equation. [START_REF] Qian | Experimental study of turbulent unconfined groundwater flow in a single fracture[END_REF] carried out laboratory experiments to study groundwater flow in a single fracture under different conditions of fracture aperture, surface roughness, and water pressure. The experimental results show that the average flow velocity (V ) could be approximated by an empirical exponential function of the hydraulic gradient (I ), which varies in the range of 0.003-0.02. The power index of the exponential function is close to 0.5. Such a V -I relationship indicates a non-Darcian turbulent flow, even though the Reynolds number is relatively low. [START_REF] Javadi | A new geometrical model for nonlinear fluid flow through rough fractures[END_REF] performed both laminar and turbulent flow simulations for a wide range of flow rates in an artificial three-dimensional fracture. They developed a new geometrical model for non-linear fluid flow through rough fractures, which suggests a polynomial expression, like the Forchheimer law, to describe the dependence of pressure drop on flow rate. Finally, this model has been evaluated with experimental results of a fracture with different geometries. A good accuracy was found between the proposed model and turbulent flow simulation results. [START_REF] Brush | Fluid flow in synthetic roughwalled fracture: Navier-Stokes, Stokes, and local cubic law simulations[END_REF] carried out simulations of fluid flow through single rough-walled fractures with various apertures using Navier-Stokes, Stokes, and local cubic law simulations. They observed that several Navier-Stokes velocity profiles have flatter peaks or noses that indicate the formation of an inertial core between the walls. They demonstrated that inertial forces can significantly influence the internal flow field within a fracture, so that the forces driving the flow field are reduced, and the overall flow rate is decreased. As the mean aperture increases, the effect of surface roughness diminishes [START_REF] Boutt | Trapping zones: The effect of fracture roughness on the directional anisotropy of fluid flow and colloid transport in a single fracture[END_REF].

In the present paper, several hydraulic tests have been carried out in a fractured rock sample, and a non-linear behavior of the flow regime has been ascertained. The study is aimed at determining the relationship between the average velocity and the hydraulic gradient for each investigated fracture-controlled flow pathway, and will serve as the first 2 Experimental setup

Characterization and preparation of block sample

The experiments have been performed on a limestone block with parallelepiped shape (0.6 × 0.4 × 0.08 m 3 ) recovered from the "Calcare di Altamura" formation, which is located in the Apulia region of southeastern Italy.

In Table 1 are shown the bulk hydraulic parameters of the limestone block. The fracture network has been made artificially through 5 kg mallet blows. The fissured system and fracture apertures on the block surfaces have been recorded with a high resolution digital camera. Subsequently, the images have been scaled and rectified using Perspective Rectifier (www.rectifiersoft.com) (Fig. 1). Profiles of discontinuities and aperture measurements have been extracted from the recorded images using edge function with the "canny" filter from the built-in Scilab Image Processing Toolbox (www.scilab.org). For each discontinuity the median profile and the aperture distribution have been determined (Table 2).

The surface of the block sample has been sealed with transparent epoxy resin [START_REF] Leven | Investigation of the effects of fractured porous media on hydraulics tests -an experimental study at laboratory scale using single well methods[END_REF]. A hole of 1 cm diameter has been opened for each discontinuity in correspondence to the boundary of the block.

The spatial position of opened ports (Fig. 2) and the illustration of construction details (Fig. 3) have been reported.

Materials and methods

In Fig. 4 the diagram of the experimental set up is shown. The sealed block sample is connected with a hydraulic circuit. Water moves from the upstream to the downstream tank and returns to the upstream tank by means of a transfer pump. A flow cell is connected to the outlet port. The sealed block and the tubes of hydraulic circuits are completely saturated. Initially, the valves "a" and "b" are closed, and the hydrostatic head in the flow cell is equal to h 0 . The ultrasonic flow velocimeter measures the snapshot flow rates that enter the sealed block. The experiment begins with the opening of the valve "a", which is reclosed when the hydraulic head in the flow cell is equal to h 1 . Finally, the hydraulic head in the flow cell is reported to h 0 through the opening of the valve "b".

The average flow rate through the sealed block can be estimated by means of the volumetric method:

Q = S 1 t 1 -t 0 (h 1 -h 0 ) (9)
where S 1 (L 2 ) represents the cross-sectional area of the flow cell. During the experiments, these values are compared with the snapshot flow rates measured by the ultrasonic velocimeter in order to check the absence of hydraulic losses due to obstructions and leaks in the hydraulic system. The experiment is repeated changing the hydraulic head h c (control head) of the upstream tank, and for each configuration of inlet-outlet ports. For different values of h c , the time, t = (t 1 -t 0 ), required to fill the flow cell from h 0 to h 1 has been registered.

The storage property of the upstream tank (S 2 ) is much higher than that of the downstream flow cell (S 2 ≫ S 1 ), therefore during the experiment the upstream hydraulic head can be considered constant. Given that the sealed block sample and the hydraulic circuit are very rigid, their compressibility can be neglected.

On the basis of these assumptions the drainage process is governed by the following equation: where h (L) is the hydraulic head of the downstream flow cell, h c (L) is the hydraulic head of the upstream tank, and

S 1 dh dt = Ŵ ( h) (h c -h) (10)
Ŵ ( h) (L 2 T -1
) is the hydraulic conductance term (Harbough, 2005) representative of both hydraulic circuit and the active fracture network configuration. Hydraulic losses within the single hydraulic circuit can be expressed according to Chezy law as:

Q = C | h| ⇒ h = 1 C 2 Q 2 (11)
where C is a characteristic coefficient related to the roughness, section and length of the tubes of the hydraulic circuit. Whereas, only for the sealed block, the h-Q relationship can be expressed through the following polynomial expression:

h = A • Q + B • Q 2 (12)
where A and B are the linear and non-linear hydraulic loss coefficients, respectively, and are related to the roughness, aperture, length and shape of fractures. Combining Eqs. ( 11) and ( 12), the conductance term representative of the whole hydraulic system assumes the following expression:

Ŵ ( h) = 2 A + A 2 + 4 B + C -2 | h| . ( 13
)
Substituting Eq. ( 13) into Eq. ( 10) and integrating the latter from t = t 0 to t = t 1 with the initial condition h = h 0 , the following equation is obtained:

S 1 -A 2 + 4(B + C -2 )(h c -h) + -A ln A 2 + 4(B + C -2 )(h c -h) -A h 1 h 0 = t 1 -t 0 . (14)
Then, fitting the experimental relationship between the time, t = (t 1 -t 0 ), and the hydraulic head of the downstream

Experimental results

Several experiments have been carried out for each in-out port configuration. The control head h c varies in the range of 0.17-1.37 m, and the average flow rates observed are in the range of 3.08 × 10 -7 -2.99 × 10 -5 m 3 s -1 . All the executed experiments show a non-linear Q-h relationship that can be well described by Eq. ( 12). Figure 5 shows the fitting method described in the previous section to estimate the linear (A) and non-linear (B) terms. The double entry Table 3 shows the estimates of A and B for each pair of ports. Figure 6 shows flow rate Q versus hydraulic head difference dh for all the experiments. If the experiments were carried out on a single fracture, a f and b f coefficients of Forchheimer Eq. ( 5) could be derived in an analytical way from Eq. ( 12). In the present case, in order to obtain Forchheimer terms for each path a numerical model has been implemented.

Starting from fracture profiles (Fig. 1), the threedimensional geometry of the fissured system has been reconstructed using the GMSH pre-processing tool (Geuzaine and Remacle, 2009). The geometry has been imported in Comsol Multiphysics ® v 4.0a (Comsol AB, 2010) using STL exchange file format (Vinciguerra and Bernabè, 2009). Furthermore, the geometry of the port holes has been modeled.

COMSOL uses the finite element scheme to solve generic partial differential equations. In particular, "Weak Form Boundary PDE Interfaces" included in "Mathematics Module" has been used to solve the following continuity equation:

∇ •   2 a f + a 2 f + 4b f ∇h ∇h   = 0. ( 15
)
Figure 7a shows the mesh of the finite element model used for numerical simulations. a f and b f are constant for the whole domain and represent the equivalent parameters for linear and non-linear hydraulic losses. The aim of the numerical model is to estimate a f and b f for each port configuration. By means of Eq. ( 9), an estimate of flow rate Q obs for steady-state conditions can be obtained in correspondence to a hydraulic head difference imposed between the inlet and outlet ports. In a similar way by means of numerical simulations, a relationship between the hydraulic head difference imposed between inlet and outlet ports and the simulated flow rate Q sim in correspondence to outlet ports can be obtained by varying a f and b f . The idea is to find the equivalent parameters a f and b f for each path that minimizes the difference between Q obs , evaluated through Eq. ( 9), and Q sim , evaluated by means of numerical model for each imposed hydraulic difference. The double entry Table 4 shows a f and b f estimated for each pair of ports. 1 -2.96e3 6.34e3 2.64e3 5.29e3 3.11e3 3.02e4 2 1.94e7 -8.04e3 3.83e3 6.30e3 3.88e3 2.99e4 3 5.77e7 1.02e8 -5.86e3 8.43e3 6.97e3 3.00e4 4 9.32e6 2.91e7 6.33e7 -6.10e3 2.97e3 2.06e4 5 5.69e7 8.34e7 1.42e8 7.15e7 -1.33e4 2.62e4 6 1.66e7 3.39e7 8.92e7 1.51e7 3.80e8 -3.91e4

Discussion

In order to analyze the experimental results, two dimensionless numbers have been evaluated: the Reynolds number and the Forchheimer number.

Reynolds number is defined as the ratio of inertial forces to viscous forces:

Re = ρ vD h µ (16)
where v represents the average velocity evaluated on the active path, and D h represents the characteristic dimension. For a fracture having small aperture with respect to its height, the characteristic dimension radius is equal to D h = w/2. Assuming that the average specific discharge is equal to q = v • w, the Reynolds number becomes:

Re = ρ 2µ q. ( 17
)
The Forchheimer number is the ratio of non-linear to linear pressure losses:

Fo = kβρ v µ . ( 18
)
Reynolds number versus Forchheimer number for all paths. According to Eq. ( 8) the Forchheimer number can be reformulated as:

Fo = b f • q a f . ( 19
)
According to [START_REF] Zeng | A Criterion for Non-Darcy Flow in Porous Media[END_REF], the Forchheimer number is recommended as a criterion for identifying non-Darcy flow because it has the advantage of clear meaning. It equals the ratio of pressure drop caused by liquid-solid interactions to that by viscous resistance, and it is directly related to the non-Darcy effect. Inertial effects dominate over viscous effects at the critical Forchheimer number (Fo > 1) [START_REF] Ruth | On the Derivation of the Forchheimer Equation by Means of the Averaging Theorme Transport in Porous Media[END_REF]. The Reynolds number is instead a dimensionless number that indicates when microscopic inertial effects become important. It is inappropriate on the macroscopic level, because microscopic inertial effects do not directly lead to macroscopic inertial effects. In fact, a high microscopic Reynolds number does not necessarily imply non-Darcian flow. Instead, Fo indicates precisely the onset of non-Darcian flow [START_REF] Ruth | On the Derivation of the Forchheimer Equation by Means of the Averaging Theorme Transport in Porous Media[END_REF]; it accounts for both velocity (v) and structure of the medium because β is structure dependent. The term β inherently contains information on the tortuosity of the flow paths that leads to changes in the microscopic inertial terms. In fact, if the structure of the medium is such that microscopic inertial effects are rare, then β will be small, and Fo will remain small until v (i.e. Re) is large. Instead, both β and Fo will be large if the structure of the medium is such that microscopic inertial effects can be expected.

In even if the paths (1-4 and 4-6) reach relatively high values of Re, they present values of Fo lower than unity.

The slope of the straight lines represents the ratio between Re and Fo that is equal to the ratio of their respective characteristic dimensions:

ζ = D h kβ . (20)
This dimensionless group ζ is characteristic of the flow path.

Relatively high values of this parameter correspond to a more linear flow behavior, because inertial effects dominate viscous ones at higher Reynolds numbers. Therefore it permits us to distinguish a different behavior of the experiments carried out with varying configurations of ports.

In particular, in Fig. 8 it is possible to distinguish a set with shallower slopes. This set has in common the outlet port 7. In Fig. 9 the shape of the fracture in correspondence to port 7 is shown. The particular shape of this fracture gives rise to a greater contact between fracture surfaces compared to the others. In fact, the path that contains this fracture presents a very high hydraulic loss.

For each path the equivalent aperture w eq has been estimated from the linear term assuming the cubic law is valid:

w eq = 3 µ ρg 12 a f . ( 21
)
This term has been compared with the average measured aperture of each path w (Fig. 9). Though w eq underestimates w, they are of the same order of magnitude.

A power law has been observed between the terms a and b (Fig. 11):

b = 0.9145 • a 2.6533 . ( 22
)
This correlation between inertial and viscous coefficients is customarily used in petroleum production engineering in order to predict high velocity well performance. Geertsma (1974) and [START_REF] Skjetne | Experiments and modeling of high-velocity pressure loss in sandstone fractures[END_REF] found similar relationships for high velocity flow in porous media and fractures.

Zimmerman and Bodvarsson (1996) analyzed flow in twodimensional rough-walled rock fractures and found a factor equal to the ratio of the cubes of w eq and w that reflects the tortuosity induced into the streamlines by the obstacles.

In the case of a single rock fracture, this factor depends only on the area of the asperity region. However, in the case study it depends on several parameters such as roughness, the areas where the rock faces are in contact with each other, fracture intersection, position and shape of the inlet and outlet ports:

τ = w 3 eq w 3 . ( 23
)
This parameter measures how much each path deviates from the parallel plate model. The pressure drops depend by the morphology of the fracture wall surfaces and on the tortuosity of the flow paths. Significant head losses may be envisaged to occur adjacent to sharp corners of a fracture where sudden change of aperture occurs [START_REF] Javadi | A new geometrical model for nonlinear fluid flow through rough fractures[END_REF].

Figures 12 and13 show τ -b and τ -a relationships, respectively. τ results correlated with b and a by means of a power function.

Conclusions

In this paper non-linear fluid flow through rock fractures was studied by means of laboratory tests and numerical modeling.

The Forchheimer equation has proven highly effective in explaining the relationship between flow rate and pressure drop, which depicts a strong inertial regime where the viscous and inertial pressure drops are controlled by v and v 2 term, respectively.

The equivalent linear and non-linear terms of Forchheimer's law have been estimated by numerical modeling.

The equivalent aperture of each flow path is determined assuming the cubic law is valid. Though it underestimates the mean measured aperture, it keeps its same order of magnitude.

A tortuosity factor τ has been determined as the ratio of the cube of the equivalent aperture and the cube of the mean measured aperture. This factor is a measure of the deviation of each flow path from the parallel plate model. In other words it measures the effects of different factors such as roughness, the contact area between fracture surfaces, fracture intersections, and the position and the form of the inlet and outlet ports.

A power law between the Forchheimer terms and τ has been detected. In complex fracture networks the tortuosity factor plays an important role in fluid flow dynamics.

The experimental results showed that the dependence of hydraulic conductivity on specific discharge cannot be neglected in fractured media. For instance, during pumping tests a linear flow model can cause errors in the determination of transmissivity in fractured rock aquifers, because much of the data collected can be non-linear due to flow occurring in transition between linear and fully turbulent flow.

On the other hand, potential errors induced by a non-linear flow model in constant pressure tests have also been recognized in the engineering literature [START_REF] Louis | Determination of in situ hydraulic parameters in jointed rock[END_REF]. [START_REF] Elsworth | Application of non-linear flow laws in determining rock fissure geometry from single borehole pumping tests[END_REF] used mathematical modeling of packer tests in fractured rock to show that calculation of transmissivity using non-Darcian constant head data can lead to underestimation errors as much as an order of magnitude.

In our study the effective hydraulic transmissivity proves to be less than 46.59 % (average value) of the Darcian (linear) flow hydraulic transmissivity. In particular, in correspondence to path 3-4 the variation reaches a maximum value equal to 59.38 %.

In pumping tests multiple pressure steps (i.e. higher flow rates resulting in a much greater differential pressure) should be used for a more accurate identification of the Darcian range and the quantification of the linear to non-linear flow relations, resulting in better transmissivity estimates as a non-linear function of the gradient.

This concept has to be taken into account in cases of anthropogenic stresses in the aquifer that might give rise to high hydraulic gradients.

Fig. 1 .

 1 Fig. 1. Block sample rectified images (0.6 × 0.4 × 0.08 m 3 ). White traces indicate the median profiles.

kFig. 2 .

 2 Fig. 2. Localization of the ports on the top surface of the block; (a) top fracture traces, (b) bottom fracture traces, and (c) ports.

Fig. 3 .

 3 Fig. 3. Illustration of a "port" construction details.

C.Figure 4 .

 4 Figure 4. Schematic diagram of the experimental setup (picture is not to scale). 516 Fig. 4. Schematic diagram of the experimental setup (picture is not to scale).

Fig. 5 .

 5 Fig. 5. Experimental results obtained for 1-2 configuration ports. (a) Control head H c versus time. (b) Average flow rate Q versus difference head evaluated as dh = H c + (h 0 + h 1 )/2. (c) Difference head versus conductance term evaluated as Eq. (13). (d) Average flow rate versus resistance term evaluated as the inverse of conductance. Dots represent the experimental values, dashed line represent the fitting of experimental values, and marked line represent the functions without the effect of circuit (C = 0).

523Figure 6 Fig. 6 .

 66 Figure 6 Flow rate Q versus hydraulic head difference dh for all the experiments. Dots represents the experimental values,

  The upper and lower triangular matrix represents the linear coefficient (a f ) and the non-linear coefficient (b f ), respectively, obtained for each pair of the ports.

  Fig. 7. (a) Finite element mesh of the numerical model. (b) A result of a simulation, color scale indicates hydraulic head, black arrows represent specific discharge.

Fig. 8 .

 8 Fig. 8. Reynolds number versus Forchheimer number for all paths.

Fig. 9 .

 9 Fig. 9. Fractured block. The fracture in correspondence to the hole 7 is shown in the lower left side.

Fig. 10 .

 10 Fig. 10. Average measured aperture (continuous line) and equivalent aperture (dashed line) for each path.

Fig. 11 .

 11 Fig. 11. Linear term a versus non-linear term b.

Fig. 12 .

 12 Fig. 12. Tortuosity factor τ versus non-linear term b.

Fig. 13 .

 13 Fig. 13. Tortuosity factor τ versus linear term a.

Table 1 .

 1 Properties of the limestone block.

	Bulk density (g cm -3 )	2.21
	Porosity (%)	0.20
	Moisture content (%)

Table 2 .

 2 Characteristic parameters of the discontinuities.

	Fracture Mean aperture Aperture range Roughness Length
	number	b m (mm)	(mm)	ε (mm)	(mm)
	1	1.496	0.811 / 2.907	0.180	119
	2	1.069	0.541 / 1.977	0.187	88
	3	1.076	0.573 / 2.784	0.187	172
	4	1.100	0.271 / 3.900	0.101	171
	5	1.040	0.811 / 2.194	0.181	141
	6	0.917	0.272 / 1.794	0.141	107
	7	0.906	0.540 / 1.510	0.010	69
	8	1.484	0.540 / 3.097	0.150	62
	9	0.703	0.540 / 3.577	0.450	155
	10	0.875	0.810 / 1.300	0.063	65
	11	1.075	0.270 / 2.901	0.152	78
	12	0.921	0.544 / 1.912	0.135	246
	13	1.130	0.541 / 1.910	0.096	258
	14	1.576	0.344 / 3.050	0.222	143

Table 3 .

 3 Double entry table. The upper and lower triangular matrix represents the linear coefficient (A) and the non-linear coefficient (B), respectively, obtained for each pair of the ports.

	1	2	3	4	5	6	7
	1 -	1.54e4	4.16e4	1.77e4	3.06e4	1.63e4	1.47e5
	2 1.53e9	-	3.71e4	2.81e4	4.11e4	2.42e4	2.06e5
	3 2.79e9	3.75e9	-	5.78e4	7.32e4	5.62e4	2.50e5
	4 7.22e8	3.00e9	3.94e9	-	3.42e4	1.45e4	1.64e5
	5 4.01e9	6.61e9	7.80e9	3.58e9	-	2.90e4	1.83e5
	6 1.11e9	2.73e9	4.22e9	9.18e8	6.71e9	-	2.28e5
	7 3.22e11 2.66e11 2.66e11 2.79e11 2.95e11 2.97e11 -

Table 4 .

 4 Double entry table.

Hydrol. Earth Syst. Sci., 16, 2511-2522, 2012 www.hydrol-earth-syst-sci.net/16/2511/2012/

www.hydrol-earth-syst-sci.net/16/2511/2012/ Hydrol. Earth Syst. Sci., 16, 2511-2522, 2012

4.97e9 3.91e9 3.74e9 2.47e9 3.53e9 5.67e9www.hydrol-earth-syst-sci.net/16/2511/2012/ Hydrol. Earth Syst. Sci., 16, 2511-2522, 2012

Acknowledgements. The authors would like to acknowledge Alberto Guadagnini, Noelle Odling and Philippe Davy for their valuable comments and suggestions that helped improve the manuscript.

Edited by: A. Guadagnini