L Cea

C Legout

F Darboux

M Esteves

G Nord

Experimental validation of a 2D overland flow model 1 using high resolution water depth and velocity data 2

Keywords: overland flow, 2D shallow water equations, bed roughness, Large Scale Particle Image Velocimetry, laser scanner

This paper presents a validation of a two-dimensional overland flow model using empirical laboratory data. Unlike previous publications in which model performance is evaluated as the ability to predict an outlet hydrograph, we use high resolution 2D water depth and velocity data to analyze to what degree the model is able to reproduce the spatial distribution of these variables. Several overland flow conditions over two impervious surfaces of the order of one square meter with different micro and macro-roughness characteristics are studied. The first surface is a simplified representation of a sinusoidal terrain with three crests and furrows, while the second one is a mould of a real agricultural seedbed terrain. We analyze four different bed friction parameterizations and we show that the performance of formulations which consider the transition between laminar, smooth turbulent and rough turbulent flow do not improve the results obtained with Manning or Keulegan formulas for rough turbulent flow. The simulations performed show that using Keulegan formula with a physically-based definition of the bed roughness coefficient, a two-dimensional shallow water model is able to reproduce satisfactorily the flow hydrodynamics. It is shown that, even if the resolution of the topography data and numerical mesh are high enough to include all the small scale features of the bed surface, the roughness coefficient must account for the macro-roughness

INTRODUCTION

Several studies on the experimental validation of overland flow models based on the two-dimensional shallow water equations have shown that, with a proper bed friction calibration and a well defined Digital Terrain Model (DTM), these kind of models are able to produce accurate predictions of the outlet hydrograph in a basin (Cea et al., 2010a,b;[START_REF] Howes | One-and two-dimensional modelling of overland flow in semiarid shrubland, jornada basin, new mexico[END_REF][START_REF] Kivva | Two-dimensional modeling of rainfall runoff and sediment transport in small catchments areas[END_REF][START_REF] Morgali | Computer analysis of overland flow[END_REF]Yan and Kahawita, 1989;Zhang and Cundy, 1989).

The evaluation of model performance in all these studies is based on the correct prediction of the outlet hydrograph generated by a rainfall event, which is of most importance for flood forecasting applications. However, a distributed model gives much more information regarding the spatial distribution of water depth and velocity, which might be used in contaminant transport and erosion models. In fact, there has been a trend in the last years to develop physically-based water erosion and solute transport models including the representation of processes at their lowest scales (Hairsine and Rose, 1992b,a;[START_REF] Jomaa | Effect of raindrop splash and transversal width on soil erosion: Laboratory flume experiments and analysis with the Hairsine ŰRose model[END_REF]Nord andEsteves, 2007, 2010;Shaw et al., 2006[START_REF] Shaw | Accounting for surface roughness in a physically-based urban wash-off model[END_REF]. These models use formulations that require a very detailed spatial characterization of the flow field (water depth, velocity and bed friction) in order to take advantage of their full potential. The problem of predicting flow transport capacity demands not only a proper sediment transport equation, but a reliable means of modelling overland flow hydraulics on rough Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M., Nord, G. (2014). Experimental validation of an 2D overland flow model using high resolution water depth and velocity data. Journal of Hydrology, 513, 142-153. DOI : 10.1016/j.jhydrol.2014.03.052

surfaces [START_REF] Abrahams | Predicting sediment transport by interril overland flow on rough surfaces[END_REF]. In this vein, several studies have been devoted

to the physical characterization of bed roughness in order to improve the prediction capability of overland flow models [START_REF] Jarvela | Flow resistance of flexible and stiff vegetation: a flume study with natural plants[END_REF][START_REF] Shaw | Accounting for surface roughness in a physically-based urban wash-off model[END_REF], and it has been recognized in previous laboratory studies involving overland flow conditions that the effective bed roughness coefficient varies substantially with the flow properties due to additional head losses generated by small scale flow features which are not resolved in the numerical models [START_REF] Fraga | Experimental study of the water depth and rainfall intensity effects on the bed roughness coefficient used in distributed urban drainage models[END_REF][START_REF] Howes | One-and two-dimensional modelling of overland flow in semiarid shrubland, jornada basin, new mexico[END_REF][START_REF] Lawrence | Hydraulic resistance in overland flow during partial and marginal surface inundation: Experimental observations and modeling[END_REF]Shaw et al., 2006;Tatard et al., 2008).

However, an accurate prediction of the outlet hydrograph does not guarantee that the hydraulics within the plot is properly computed, and so such an approach is not a thorough model validation. One of the main conclusions of the work of Tatard et al. (2008) was that the Froude number was not predicted by any of the three models tested on a 40 m 2 plot and even the general pattern was missed.

This highlights the need to evaluate the performance of two-dimensional overland flow models in terms of their capability to reproduce the spatial distribution of water depth and velocity at the same time. Until now, very few studies have made such evaluation. Tatard et al. (2008) compared measured and modelled depths and velocities at selected points without performing a comparison on the whole plot. More recently, [START_REF] Mügler | Comparison of roughness models to simulate overland flow and tracer transport experiments under simulated rainfall at plot scale[END_REF] intended to make an evaluation of what occurs within the plot by comparing measured and modelled tracer breakthrough curves of point-source injection of tracer under rainfall conditions. Such an approach aims at assessing if the actual travel time of a tracer is well reproduced by the model. However, it is a spatially integrated approach that does not allow the direct comparison of observed and modelled depth-velocity maps at the mesh size scale. To the best of our knowledge, the work presented in this paper provides

Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M., Nord, G. (2014)

∂q x ∂t + ∂ ∂x q 2 x h + gh 2 2 + ∂ ∂y q x q y h = -gh ∂z b ∂x - τ b,x ρ (2) ∂q y ∂t + ∂ ∂x q x q y h + ∂ ∂y q 2 y h + gh 2 2 = -gh ∂z b ∂y - τ b,y ρ (3)
where z b is the bed elevation, h is the water depth, (q x , q y ) are the two components The bed friction stress in Equations (2-3) can be expressed as:

τ b,x ρ = C f |U|U x τ b,y ρ = C f |U|U y (4)
where |U| is the modulus of the depth-averaged horizontal velocity, (U experiments, which will be referred to as sinusoidal mould (Figure 1) and agricul- mental conditions the water depth was less than 10 mm in the whole mould. The Reynolds number in the main furrows varies between 1000 and 12000, depending on the test case and on the furrow considered, while the Froude number varies between zero and 5, although in most of the mould it is larger than one.

The agricultural mould is a square of 0.60 m by 0.60 m representative of a seedbed. The mould is a perfect reproduction of a natural seedbed obtained with a methodology similar to the one presented in [START_REF] Kamphorst | Validation of a numerical method to quantify depression storage by direct measurements on moulded surfaces[END_REF]. All the experiments were performed on that single reproduction. It exhibits two main parallel furrows oriented along the X-direction. The width of these furrows is of the order of 0.10-0.15 m (Figure 2e). The surface of the agricultural mould is made of plaster and is therefore much smoother than that of the sinusoidal mould. On the other hand, in the agricultural mould there is a macro-roughness formed by small changes in surface topography, as a result of different sized soil aggregates present in the actual seedbed. In order to determine how the representation of these small supplied on the whole width of the upper part of the surface. However, the water systematically enters the mould through the two main furrows, as indicated by the arrows in Figure 2a. Two combinations of flow rates and slopes in the X-direction on the agricultural mould have been used for model validation (Table 2). These

were selected in such a way that water flows either in two separated furrows or in two connected furrows. In this case, only the velocity field was available for model validation.

Datasets

For evaluation and comparison purposes, all measurements were performed at steady state. The flow rates were measured separately at the outlet of each furrow by automatically weighing the volume of water for 20 s to 2 min depending on the flow rate for the sinusoidal mould. For the agricultural mould, the total flow rates were measured by manual sample collection at the outlet every 2 min.

The maps of flow depth in the sinusoidal mould were obtained after subtraction of two elevation maps acquired successively (i.e., subtraction of the water surface map from the DTM without water). Both elevation datasets were measured using the instantaneous profile laser scanner described in [START_REF] Darboux | An instantaneous-profile laser scanner to measure soil surface microtopography[END_REF] according to the protocol presented in [START_REF] Legout | High spatial resolution mapping of surface velocities and depths for shallow overland flow[END_REF]. The spatial resolution was 0.5 mm in X-and Y-directions. (agricultural mould). Table 3 shows the properties of the different grids used in the mesh convergence analysis.

The topography of the sinusoidal mould is quite smooth and regular and therefore, its main features are well captured even with a rather coarse mesh. On the other hand the agricultural mould has a more irregular topography, including soil aggregates up to 10-20 mm in size. In this case the mesh size has an important impact on the accuracy of the bed elevation definition. The spatial resolution of meshes A1M2, A1M3 and A1M4 (Table 3) is lower than that of the original DTM and therefore, there is a loss of accuracy in the bed topography used in the numerical model.

In both the sinusoidal and the agricultural terrains, the results obtained with the two finest meshes (S1M1 and S1M2 in the sinusoidal mould and A1M1 and A1M2 in the agricultural mould) are almost identical, and show similar flow patterns to those obtained with meshes S1M3 and A1M3 respectively (Figures 3 and4). Nonetheless, the maximum water depths and velocities are slightly diffused with the latter meshes. The results obtained with the meshes S1M4 and A1M4 are clearly very diffusive and are not able to capture the maximum values of the water depth and velocity, although the general flow pattern can still be identified.

Therefore, we have decided to work with the mesh S1M2 in the sinusoidal mould and A1M2 in the agricultural mould, which have a spatial resolution of 5 mm and 4 mm respectively. These meshes show a good compromise between computation time and accuracy on model output (Table 3). Results obtained with the finest meshes (S1M1 and A1M1) do not show significant differences while the computation time increases considerably.

Comment citer ce document : Cea, L., Legout

Keulegan formulation applied to the sinusoidal mould

The surface of the sinusoidal mould is made of sand grains with a diameter within 0.5 and 1.0 mm. In non-cohesive beds without bedforms, river hydraulics manuals recommend to use a roughness height within 1 and 3 times the diameter of the sediment [START_REF] Garcia | Sedimentation Engineering: Processes, Measurements, Modeling, and Practice[END_REF]. Based on this approximation, we have used a roughness height of k s = 1.5 mm in the simulations of the sinusoidal mould test cases. It should be noticed that this approximation is recommended for river flows and not specifically for overland flows, where the water depth in some regions might be of the same order of magnitude as the roughness height. Nonetheless, as it will be shown in the following this bed roughness coefficient gives very satisfactory results in the four sinusoidal mould test cases.

Using Keulegan formulation with k s = 1.5 mm the overall experimentalnumerical agreement is very satisfactory . The model predicts correctly the number of furrows which contribute to drain the inlet discharge (3 furrows in test cases S1 and S2, two furrows in test case S3, and one single furrow in test case S4). The prediction of the water depth field is very good in the four test cases. The root mean square errors on the water depth are 0.6 mm, 0.7 mm, 0.7 mm and 0.9 mm for the test cases S1 to S4 respectively (Table 4), which are of the same order of magnitude as the diameter of the sediment which covers the mould surface. In all test cases the error on the water depth is less than 1 mm in more than 90% of the measurement points, while it is greater than 2 mm in less than 1% of the points (Figure 9). The accuracy on the prediction of the hydraulic jump which appears in the test cases S1 and S4 at the location y = 0.78 m in the profile x = 0.8 m is remarkable (Figure 11). certainty for two reasons. Firstly, the LSPIV technique gives abnormally small velocities in regions where the water depth is very small, because the tracers occasionally rest attached to the rough bed surface. This occurs typically at the flow inlet and in the crests between furrows. In order to compare the velocities only the points in which the unit discharge is larger than 0.5 cm 2 /s were considered.

This value was adopted because it defines rather precisely the shape of the furrows (Figures 5678), and it excludes the regions where the LSPIV velocity has large measurement errors. Another source of uncertainty in the comparison of velocity data is the fact that LSPIV measures the surface velocity, while the numerical model computes the depth-averaged velocity. In order to estimate the surface velocity from the model output, the ratios mentioned in [START_REF] Legout | High spatial resolution mapping of surface velocities and depths for shallow overland flow[END_REF] and in section 3.2 were adopted. With these criteria, the root mean square errors on the velocity are 0.089 m/s, 0.058 m/s, 0.083 m/s and 0.065 m/s respectively for the test cases S1 to S4 (Figure 10).

It should be noted that in the test case S1 (Figure 5), there is a clear disagreement between numerical predictions and experimental data in the third furrow, which is the one with the smallest water depth. The experimental measurements of water depth and velocity reveal a region between x = 0.1 m and x = 0.4 m where the flow is locally accelerated (higher velocity and lower water depth). The same effect appears in the test case S3 (Figure 7), but in this case in the second furrow (which is in this case the one which presents the smallest water depths).

This local acceleration does not appear in the numerical model, as it can be observed in the velocity fields shown in Figures 5 and8, and it is the reason why the root mean square error is larger in these test cases than in the cases S2 and the LSPIV velocity is clearly higher than the numerical velocity. We have tried unsuccessfully to improve the velocity predictions in these regions using different bed friction formulations (see section 4.3). We have not found a clear physical explanation for this effect.

Other bed friction formulations applied to the sinusoidal mould

The results shown in the previous section indicate that Keulegan formulation with a definition of the roughness height based on the physical properties of the mould surface (in our case the sediment diameter) performs satisfactorily in the sinusoidal mould. As roughness formulations are known to affect significantly the model performance [START_REF] Mügler | Comparison of roughness models to simulate overland flow and tracer transport experiments under simulated rainfall at plot scale[END_REF], we have explored the possibility of improving the predictions of the model by using other bed friction parameterizations.

For that purpose we have considered the standard bed friction formulations shown in Table 1. Three of these formulations (Keulegan, Manning and Darcy with a constant friction factor) were developed for rough turbulent flow. The forth formulation is based on the formula of Colebrook-White, and considers the transitions between laminar flow, smooth turbulent flow and rough turbulent flow. In this latter formulation the laminar bed friction formula is used if the local Reynolds number is smaller than 2000 (4U h/ν < 2000), otherwise the Colebrook-White equation is applied. In all cases the bed roughness parameter was calibrated for the test case S1 and then used in the test cases S2, S3 and S4.

After calibration of their respective friction parameters, all the parameterizations give a similar agreement with the experimental data (Tables 4 and4).

Their performance considering the water depth and velocity spatial distributions is rather satisfactory, specially regarding the water depth. The agreement in the velocity field is slightly better when using Keulegan or Manning formulations,

Agricultural mould

The most significant difference between the sinusoidal and the agricultural moulds is that the latter one includes soil aggregates which create a certain level of macro-roughness. As shown in Figure 2, the size and location of the soil aggregates is well captured by the numerical meshes with a spatial resolution of 2 mm and 4 mm (A1M1 and A1M2), but they are reduced as the grid resolution is reduced. The macro-roughness generated by these soil aggregates induces small scale accelerations and decelerations in the flow field, causing additional head losses. One of the purposes of the agricultural terrain test cases is to analyze if a high resolution shallow water model with a micro-roughness coefficient defined from the bed material (in this case plaster) is able to reproduce the experimental velocity field.

Test case A1 has been used for model calibration and to analyze how the bed friction coefficient should be modified in order to account for the macroroughness. Test case A2 has been used for model validation, using the same roughness parameterization as in case A1.

Figure 12 shows that the velocity field computed using a bed roughness coefficient for plaster (either k s = 1 mm in Keulegan formulation or n = 0.013 s m -1/3 in Manning formulation) are far from the experimental ones, even if the computations were done with a high resolution mesh which includes all the small scale details of the topography. The computed velocities are too large compared with the LSPIV measurements. The root mean square error obtained with Keulegan and Manning formulas using these parameters are respectively 0.095 m/s and 0.100 m/s (Table 5).

Model outputs is improved significantly if a larger bed friction coefficient is used. After calibration, the coefficients which produce the best numericalexperimental agreement were found to be k s = 30 mm and n = 0.050 s m -1/3 , the root mean square error for the velocity being 0.040 m/s for both cases. The velocity fields computed with these parameters show a good global agreement with the LSPIV data (Figure 12). Both, experimental and numerical results show a shallow and narrow arm which develops to the left of the upper furrow and joins the lower furrow further downstream.

It can be noticed in Figure 12 shown in Figure 13.

The same roughness parameterizations applied to the test case A2 gives satisfactory results (Figure 13). Although the quantitative agreement between numerical and experimental data is not as good as in test case A1 (Table 5), the differences are of the same order as the uncertainty on the experimental surface velocity (section 3.2). This increase in the root mean square error is to some extent explained by the gaps in the experimental velocity data which can be observed in Figure 13. Due to the extremely low water depths and the high surface roughness of the agricultural mould, there are some isolated spots in which tracers could not move and thus velocity was not correctly measured. This introduces spurious disagreements with the numerical model which are reflected in the root mean square error. Nonetheless, visual comparison of the velocity fields in Figure 13 shows a similar global agreement to that of test case A1.

Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M., Nord, G. (2014). Experimental validation of an 2D overland flow model using high resolution water depth and velocity data. Journal of Hydrology, 513, 142-153. DOI : 10.1016/j.jhydrol.2014.03.052

CONCLUSIONS

The performance of an overland flow model based on the two-dimensional shallow water equations has been analyzed using high-resolution distributed measurements of water depth and velocity. Several flow conditions over two impervious moulds with different micro and macro-roughness characteristics have been used for model validation.

The flow in the sinusoidal mould, which represents a terrain with a regular topography and no macro-roughness, is properly represented by the two-dimensional shallow water equations with a constant bed friction coefficient. This latter can be related to the micro-roughness characteristics of the bed surface. Manning and Keulegan formulations produce results with similar accuracy. More complex formulations which consider the transitions between laminar, smooth turbulent and rough turbulent flow conditions, do not improve model performance.

In terrains with a certain level of macro-roughness, as it is the case of the agricultural mould with soil aggregates, the two-dimensional shallow water equations with a bed friction coefficient defined from the micro-roughness of the bed surface are not able to accurately reproduce the flow field. This is because the soil aggregates which define the macro-roughness produce small scale accelerations and decelerations in the flow field, which generate additional head losses. Even if the mesh size is smaller than the soil aggregates, these head losses cannot be captured by the model equations. In order to correctly reproduce the flow hydrodynamic, the bed roughness coefficient must account for the macro-roughness of the terrain generated by the presence of soil aggregates.

In surface (l ∆ ) by the simple relation (k s = l ∆ ). In the test cases analyzed in this work a proper definition of l ∆ is the maximum residual obtained after filtering the original DTM with a filter size of 10 mm, and it is related to the size of the soil aggregates. The problem with this approach when applied to larger scales is that a very fine spatial resolution is required in order to apply a 10 mm filter to the DTM, which limits the applicability of a macro-roughness based definition of k s .

These results show that using a physically-based definition of the bed roughness coefficient, a two-dimensional shallow water model is able to reproduce the velocity and water depth patterns under runoff conditions. In the test cases analyzed in this paper the average errors on the water depth and velocity were 0.7 mm and 0.07 m/s respectively. This justifies the use of the shallow water equations in physically-based erosion and solute transport models that require a very detailed spatial characterization of the flow field in order to take advantage of their full potential. Regarding the applicability of the model to larger scales, and thus more relevant for hillslope hydrology understanding, the limiting factor is the DEM and mesh resolution required for the computations. A spatial resolution of the order of 10 mm is required to capture the flow patterns and the range of depth and velocity values. Applying such a high resolution shallow water model to plot scales of a few hundreds of square meters would lead to reasonable computation times (of the order of 5 h in a standard laptop for a 100 m 2 plot). Such resolution is also consistent with what can be obtained from stereophotogrammetry or laser scanner techniques in field conditions.

Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M., Nord, G. (2014)

Formulation

Friction factor parameterization Darcy (f ) 4). Axes units in meters.

C f = f 8 Keulegan (k s ) C f = 0.16 (ln 11h/k s ) -2 Manning (n) C f = g n 2 h 1/3 Colebrook-White (k s ) C f = 0.
Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M., Nord, G. (2014) 4). Axes units in meters.

Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M., Nord, G. (2014) 4). Axes units in meters.

64

  of the unit discharge, (τ b,x , τ b,y ) are the two horizontal components of the bed fric-65 tion stress, ρ is the water density and g is the gravity acceleration. Turbulent shear 66 stresses have been neglected in Equations (2-3) because under runoff conditions, 67 with water depths of a few millimetres, their influence in the momentum equations 68 is negligible compared to bed friction. 69 Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M., Nord, G. (2014). Experimental validation of an 2D overland flow model using high resolution water depth and velocity data. Journal of Hydrology, 513, 142-153. DOI : 10.1016/j.jhydrol.2014.03.052

102

  tural mould (Figure2a). The sinusoidal mould is a square of 1.20 m by 1.20 m 103 which consists of three parallel crests and furrows with a sinusoidal shape, with 104 a wave length of 0.20 m and an amplitude of 0.01 m. The DTM was measured 105 with a spatial resolution of 0.5 mm. The surface of the sinusoidal mould is made 106 of a thin layer of sand with a mean diameter of approximately 0.5-1.0 mm, which 107 is glued to an impervious surface. Therefore, there is no infiltration in the experi-108 ments. We will refer by micro-roughness to the bed roughness linked to the sand 109 texture. The water supply was located on the upper part of the surface. The inlet 110 was 0.10 m wide and was set perpendicular to the direction of the largest slope 111 on each test case. Four different combinations of flow rates and slopes on the si-112 nusoidal mould have been considered for model validation (

  topographic changes is affected by the DTM resolution, we have filtered the original DTM using average filters of different size, and we have defined the residual length at each point of the mould as the difference between the original and filtered DTM's. The spatial resolution of the original agricultural mould DTM is 2 mm. A filter of 4 mm produces almost no effect on the original DTM, which confirms that this resolution is high enough to capture all the topography features. A filter of 8 mm gives an average residual length of 0.35 mm and a maximum residual length of 6 mm, which is of the order of magnitude of the water depth. If the filter size is increased to 16 mm the maximum residual length increases to 12 mm. Larger filters introduce a smoothing effect in the macrotopography and distort the definition of the two main longitudinal furrows (Figure2c and d). The water was Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M., Nord, G. (2014). Experimental validation of an 2D overland flow model using high resolution water depth and velocity data. Journal of Hydrology, 513, 142-153. DOI : 10.1016/j.jhydrol.2014.03.052

  The uncertainty was experimentally assessed to 0.5 mm with maximum values of 1 mm. It should be remarked that the maps of flow depths were not measured on the agricultural mould. The maps of surface velocity measurements were obtained by Large Scale Particle Image Velocimetry. Based on the acquisition of successive images of tracers flowing on the surface, the LSPIV technique estimates the movement of tracer particles between image pairs using a cross-correlation analysis. More details on Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M., Nord, G. (2014). Experimental validation of an 2D overland flow model using high resolution water depth and velocity data. Journal of Hydrology, 513, 142-153. DOI : 10.1016/j.jhydrol.2014.03.052 the adaptation of this well-known technique to the case of very shallow flows such 168 as overland flow are presented in Legout et al. (2012). The spatial resolution of 169 the velocity map was 5 mm. The ratios between surface velocities obtained with 170 LSPIV and the average velocity on the water column were estimated in Legout 171 et al. (2012) as 0.8 and 0.9 respectively on the sinusoidal and agricultural mould. 172 In the following section, the velocity measurements were corrected by those ratios 173 for the comparison between the experimental surface velocity and the modelled 174 depth-averaged velocity. The uncertainty of the surface velocities was experimen-175 tally assessed to 0.03 m/s with maximum values of roughly 0.10 m/s.

  we analyze the performance of the overland flow model pre-178 sented in section 2 against the experimental data described in section 3. Model 179 performance is evaluated by comparing the spatial distribution of water depth 180 and velocity. In order to quantify the numerical-experimental agreement, the root 181 mean square error is computed for each simulation. As mentioned in section 3, 182 in the agricultural mould there are no available water depth measurements and 183 therefore, only the velocity field is used for comparison.

  is a critical parameter in any numerical model, as it has im-186 portant consequences on the accuracy of the results, as well as on the compu-187 tation time required to solve the model equations. Considering that the level of 188 macro-roughness of the two terrains considered in this study is very different, an 189 independent mesh convergence analysis was performed for each terrain, using for 190 that purpose the flow conditions of the test cases S1 (sinusoidal mould) and A1 191 Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M., Nord, G. (2014). Experimental validation of an 2D overland flow model using high resolution water depth and velocity data. Journal of Hydrology, 513, 142-153. DOI : 10.1016/j.jhydrol.2014.03.052

  The comparison of the measured and computed velocities has a larger un-Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M., Nord, G. (2014). Experimental validation of an 2D overland flow model using high resolution water depth and velocity data. Journal of Hydrology, 513, 142-153. DOI : 10.1016/j.jhydrol.2014.03.052

S4.

  This discrepancy is represented in Figure 10 by the cloud of points in which Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M., Nord, G. (2014). Experimental validation of an 2D overland flow model using high resolution water depth and velocity data. Journal of Hydrology, 513, 142-153. DOI : 10.1016/j.jhydrol.2014.03.052

  Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M., Nord, G. (2014). Experimental validation of an 2D overland flow model using high resolution water depth and velocity data. Journal of Hydrology, 513, 142-153. DOI : 10.1016/j.jhydrol.2014.03.052 although the performance of Darcy and Colebrook-White formulations is still acceptable. The main advantage of Keulegan parameterization is that the roughness height is more easily related to the roughness characteristics of the terrain than the Manning or Darcy coefficients.

  that the numerical model gives some locally high velocity values which do not appear in the experimental data. These local velocity peaks are due to the highly irregular topography created by the soil aggregates. The energy losses induced by these high frequency accelerations are not well represented by the shallow water equations with the standard bed friction formulations, even if a very detailed grid and DTM are used in the computations. This explains the high value of the bed friction coefficient found in the calibration (k s = 30 mm and n = 0.050 s m -1/3 ).The roughness coefficient in Keulegan formulation might be related to the macro-roughness of the bed surface. For this purpose, the roughness height in Keulegan formulation was defined on a physical basis as k s = Al ∆ , where l ∆ is a local measure of the macro-roughness length scale, and A is a constant of order 1. The level of macro-roughness in the agricultural mould is not homogeneous, the spatial density and size of the soil aggregates in the upper furrow being lower Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M., Nord, G. (2014). Experimental validation of an 2D overland flow model using high resolution water depth and velocity data. Journal of Hydrology, 513, 142-153. DOI : 10.1016/j.jhydrol.2014.03.052 than in the lower furrow. Therefore, a spatially variable macro-roughness length scale l ∆ was defined from the residuals (∆z) obtained after filtering the original DTM with an average filter of size ∆. After trying several definitions of l ∆ and filter sizes, a suitable definition for the roughness height was obtained with A = 1, and l ∆ = max(∆z), where the residuals ∆z are computed using an average filter of size ∆ = 10 mm. The velocity results obtained using this macro-roughness based definition for the roughness height (k s = max(∆z)) are shown in Figure 12. There is not a significant improvement relative to the results obtained with a constant roughness coefficient (k s = 30 mm), but in this case the values of k s vary locally according to the physical characteristics of the macro-roughness, as

  the agricultural mould the bed roughness height used in Keulegan formulation (k s ) could be related to the local macro-roughness length scale of the bed Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M., Nord, G. (2014). Experimental validation of an 2D overland flow model using high resolution water depth and velocity data. Journal of Hydrology, 513, 142-153. DOI : 10.1016/j.jhydrol.2014.03.052

Figure 1 :

 1 Figure 1: Representation of the sinusoidal mould.

Figure 2 :

 2 Figure 2: Representation of the agricultural seedbed mould after filtering the original DTM with several filter sizes (∆). (a) Original DTM, (b) ∆ = 4 mm, (c) ∆ = 8 mm (d) ∆ = 16 mm , and (e) cross section at x = 0.3 m computed from the filtered DTM's.

Figure 3 :

 3 Figure 3: Mesh convergence in the sinusoidal mould. Water depth (m) computed with the meshes S1M1 (top-left), S1M2 (top-right), S1M3 (bottom-left) and S1M4 (bottom-right). Axes units in meters.

Figure 4 :Figure 5 :

 45 Figure 4: Mesh convergence in the agricultural mould. Black arrows indicate inflow locations. Velocity (m/s) computed with the meshes A1M1 (top-left), A1M2 (top-right), A1M3 (bottomleft) and A1M4 (bottom-right). Axes units in meters.

Figure 6 :

 6 Figure 6: Experimental (left) and numerical (right) water depth (m) and surface velocity (m/s) fields for the test case S2. Numerical results were obtained with Keulegan formulation and k s = 1.5 mm. The black contour defines the region used to compute the root mean square error for the velocity (Table4). Axes units in meters.

Figure 7 :

 7 Figure 7: Experimental (left) and numerical (right) water depth (m) and surface velocity (m/s) fields for the test case S3. Numerical results were obtained with Keulegan formulation and k s = 1.5 mm. The black contour defines the region used to compute the root mean square error for the velocity (Table4). Axes units in meters.

Figure 8 :

 8 Figure 8: Experimental (left) and numerical (right) water depth (m) and surface velocity (m/s) fields for the test case S4. Numerical results were obtained with Keulegan formulation and k s = 1.5 mm. The black contour defines the region used to compute the root mean square error for the velocity (Table4). Axes units in meters.

Figure 9 :

 9 Figure 9: Computed vs. modelled water depths in the sinusoidal mould. Test cases S1 (topleft), S2 (top-right), S3 (bottom-left) and S4 (bottom-right). Lines above and below the 1:1 line correspond to the maximum experimental uncertainty of 1 mm.

Figure 10 :

 10 Figure 10: Computed vs. modelled velocities in the sinusoidal mould. Test cases S1 (top-left), S2 (top-right), S3 (bottom-left) and S4 (bottom-right). Lines above and below the 1:1 line correspond to the maximum experimental uncertainty of 0.10 m/s.

Figure 12 :

 12 Figure 12: Velocity field (m/s) in the agricultural mould in the test case A1 using different bed roughness parameterizations. Black arrows indicate inflow locations. Experimental velocity (topleft) vs. numerical velocity computed with k s = max ∆z (top-right), k s = 1 mm (middle-left), k s = 30 mm (middle-right), n = 0.013 s m -1/3 (bottom-left) and n = 0.050 s m -1/3 (bottomright).

Figure 13 :

 13 Figure 13: Velocity (m/s) field in the agricultural mould in the test case A2 using different bed roughness parameterizations. Black arrows indicate inflow locations. Experimental velocity (topleft) vs. numerical velocity computed with k s = max ∆ z (top-right) and k s = 30 mm (bottomleft). The bottom-right figure shows the spatial distribution of the roughness height computed from k s = max ∆z.

Table 2

 2

	). The slopes

Table 1 :

 1 . Experimental validation of an 2D overland flow model using high resolution water depth and velocity data. Journal of Hydrology, 513, 142-153. DOI : 10.1016/j.jhydrol.2014.03.052 Bed friction formulations considered in the numerical simulations. f : Darcy's friction factor; k s : roughness height; n: Manning coefficient; ν: kinematic viscosity of water (10 -6 m 2 /s).

Table 2 :

 2 Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M.,Nord, G. (2014). Experimental validation of an 2D overland flow model using high resolution water depth and velocity data. Journal of Hydrology, 513, 142-153. DOI : 10.1016/j.jhydrol.2014.03.052 Summary of experimental test cases used for model validation.

			16 -ln	k s 11h	+	11.9U h C f 2.5ν	-2
	Laminar	C f =	3ν |U|h
	Shaw, S.B., Walter, M.T., Steenhuis, T.S., 2006. A physical model of particulate
	wash-off from rough impervious surfaces. Journal of Hydrology 327, 618-626.
	Tatard, L., Planchon, O., Wainwright, J., Nord, G., Favis-Mortlock, D., Silvera,
	N., Ribolzi, O., Esteves, M., Huang, C.H., 2008. Measurement and modelling
	of high-resolution flow-velocity data under simulated rainfall on a low-slope
	sandy soil. Journal of Hydrology 348, 1 -12.
	Toro, E.F., 2001. Shock-capturing Methods for Free-Surface Shallow Flows. Wi-
	ley, Chichester, West Sussex PO19 1UD, England.
	Yan, M., Kahawita, R., 1989. Simulating the evolution of non-point source pollu-
	tants in a shallow water environment. Chemosphere 67, 879-885.
	Zhang, Cundy, 1989. Modeling of two-dimensional overland flow. Water Re-
	sources Research 25, 2019-2035.

Table 3 :

 3 Mesh convergence analysis for the sinusoidal and agricultural moulds. The results obtained with meshes S1M1 and A1M1 are used as a reference to compute the root mean square error in water depth (σ h ) and velocity (σ v ) on the coarser meshes. ∆x is the mesh size, N nod is the number of elements in the mesh, and T c is the computation time with an Intel Core i7 -1.60Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M.,Nord, G. (2014). Experimental validation of an 2D overland flow model using high resolution water depth and velocity data. Journal of

	GHz.
	Test case Mesh	∆x	N nod	T c	σ h	σ v
	S1	S1M1 2 mm 170817 5400 s	-	-
	S1	S1M2 5 mm	38971	340 s 0.17 mm 0.018 m/s
	S1	S1M3 10 mm	9748	35 s	0.28 mm 0.034 m/s
	S1	S1M4 20 mm	2467	4 s	0.67 mm 0.065 m/s
	A1	A1M1 2 mm 115840 4200 s	-	-
	A1	A1M2 4 mm	28960	600 s	-	0.013 m/s
	A1	A1M3 8 mm	6400	60 s	-	0.031 m/s
	A1	A1M4 16 mm	1600	8 s	-	0.265 m/s

Hydrology, 513, 142-153. DOI : 10.1016/j.jhydrol.2014.03.052

Table 4 :

 4 Performance of different bed friction parameterizations in the test cases S1, S2, S3 and S4, evaluated as the root mean square error of the water depth (σ h ) and velocity (σ v ).

	Parameterization	σ h (mm)		σ v (m/s)
		S1 S2 S3 S4	S1	S2	S3	S4
	Keulegan (k s = 1.5 mm)	0.6 0.7 0.7 0.9 0.089 0.058 0.083 0.065
	Manning (n = 0.015 s m -1/3 )	0.6 0.7 0.6 0.9 0.094 0.070 0.084 0.068
	Darcy (f = 0.16)	0.6 0.8 0.8 1.1 0.099 0.063 0.088 0.100
	Colebrook-White (k s = 1.5 mm) 0.8 0.7 0.6 1.0 0.098 0.067 0.090 0.110

Table 5 :

 5 Performance of different bed friction parameterizations in the test cases A1 and A2 in the agricultural mould, evaluated as the root mean square error of the velocity.

	Parameterization Friction factor	A1	A2
	Keulegan	k s = 1 mm	0.095 m/s 0.140 m/s
	Keulegan	k s = 30 mm	0.040 m/s 0.062 m/s
	Manning	n = 0.013 s m -1/3 0.100 m/s 0.150 m/s
	Manning	n = 0.050 s m -1/3 0.040 m/s 0.071 m/s
	Keulegan	k s = max ∆ z	0.039 m/s 0.062 m/s

Comment citer ce document : Cea, L., Legout, C., Darboux, F., Esteves, M., Nord, G. (2014). Experimental validation of an 2D overland flow model using high resolution water depth and velocity data. Journal of Hydrology, 513, 142-153. DOI : 10.1016/j.jhydrol.2014.03.052


H Thorsdottir

I M Verloop
email: verloop@irit.fr

M Verloop

Markov-modulated M/G/1-type queue in heavy traffic and its application to time-sharing disciplines

Keywords: Markov-modulation, Heavy traffic, Discriminatory processor sharing, Single-server queue

This paper deals with a single-server queue with modulated arrivais, ser vice requirements and service capacity. In our first result, we derive the mean of the total workload assuming generally distributed service requirements and any service discipline which does not depend on the modulating environment. We then show that the workload is exponentially distributed under heavy-traffic scaling. In our second result, we focus on the discriminatory processor sharing (DPS) discipline. Assurning exponential, class-dependent service requirements, we show that the joint distribution of the queue lengths of different customer classes under DPS undergoes a state-space collapse when subject to heavy-traffic scaling. That is, the limiting distribution of the queue-length vector is shown to be exponential, tirnes a deterrninistic vector. The distribution of the scaled workload, as derived for general service disciplines, is a key quantity in the proof of the state-space collapse.

Introduction

Markov-modulation is a way to formalize the embedding of queues in a random environment. The parameters of the queue in question, typically arrival rates, service requirements or both, are governed by an external Markov chain, thereby creating an extra layer of randomness around the stochastic queueing process. For classical results on Markov-modulated single-server queues with the first-come-first-serve (FCFS) discipline see, for example, [START_REF] Neuts | Matrix-geometric Solutions in Stochastic Models, an Algorithmic Approach[END_REF][START_REF] Purdue | The M/M/1 queue in a Markovian environment[END_REF][START_REF] Regterschot | The queue M/G/1 with Markov-modulated arrivals and service[END_REF]. Recent work on systems in a Markov-modulated environment can be found in, for example, [START_REF] Fiems | Markov-modulated stochastic recursive equations with applications to delaytolerant networks[END_REF][START_REF] Kim | A single server queue with Markov modulated service rates and impatient customers[END_REF][START_REF] Zorine | On ergodicity conditions in a polling model with Markov modulated input and statedependent routing[END_REF].

In this paper we will analyse a modulated queue under a heavy-traffic scaling, that is, evaluate the system at its critical load. It is a well-known result from the literature on single-server queues that, under fairly general conditions [START_REF] Kingman | The single server queue in heavy traffic[END_REF], the steadystate distributions of the appropriately scaled queue length and workload become exponential when the critical load is approached. This property has been seen to carry over to certain systems where arrivals and service times are modulated by an external Markov process; see [START_REF] Asmussen | The heavy traffic limit of a class of Markovian queueing models[END_REF][START_REF] Dimitrov | Single-server queueing system with Markov-modulated arrivals and service times[END_REF][START_REF] Falin | Heavy traffic analysis of M/G/1 type queueing systems with Markov-modulated arrivals[END_REF]. In fact, [START_REF] Asmussen | The heavy traffic limit of a class of Markovian queueing models[END_REF] establishes an even stronger result: convergence of the queue-length process to a reflected Brownian motion. Multi-class single-server queues under a heavy-traffic scaling have been studied in, for example, [START_REF] Katsuda | State-space collapse in stationarity and its application to a multiclass single-server queue in heavy traffic[END_REF] for FCFS with feedback routing, [START_REF] Izagirre | Heavy-traffic analysis of a multi-class queue with relative priorities[END_REF] for the discriminatory random order of service discipline and [START_REF] Grishechkin | On a relationship between processor sharing queues and Crump-Mode-Jagers branching processes[END_REF][START_REF] Verloop | Heavy-traffic analysis of a multiple-phase network with discriminatory processor sharing[END_REF] for the discriminatory processor sharing (DPS) policy. In particular, [START_REF] Izagirre | Heavy-traffic analysis of a multi-class queue with relative priorities[END_REF][START_REF] Verloop | Heavy-traffic analysis of a multiple-phase network with discriminatory processor sharing[END_REF] show that the steady-state queue-length vector undergoes a socalled state-space collapse and converges to an exponentially distributed variable, times a deterministic vector. The cited multi-class results under heavy-traffic scaling are all for non-modulated systems. In light of this, we will in this paper put special emphasis on a modulated multi-class single-server queue, and the limiting steady-state queue-length distribution is derived.

While there is little ambiguity in how arrival rates are modulated, there are in the literature typically two ways in which to modulate the service rates. One can (i) let the departure rate be continuously modulated throughout a customer's service, and the other approach is to (ii) let a customer's service requirement distribution be based on the state of the environment when it arrives and remain the same until the customer departs. We note that by adapting the number of different customer classes, the fixed service requirements of case (ii) can be seen as a special case of the continuously modulated requirements (i); we further elaborate on this later in the paper in Remark 1 in Sect. 2.

The way the load or traffic intensity for modulated queues is defined goes hand in hand with the way the service rates are affected by the environment. In case (i), the load is typically the average of arrival rates (where the averaging is with respect to the equilibrium distribution of the environment), say λ ∞ , divided by the average of service rates, say μ ∞ (see, for example, [START_REF] Neuts | Matrix-geometric Solutions in Stochastic Models, an Algorithmic Approach[END_REF]). In case (ii), the load is taken as the average over the arrival rate times the mean service requirement, say, λ d /μ d per state d of the environment (see for example, [START_REF] Dimitrov | Single-server queueing system with Markov-modulated arrivals and service times[END_REF][START_REF] Regterschot | The queue M/G/1 with Markov-modulated arrivals and service[END_REF]). The two load definitions represent different scenarios. In particular, when load (ii) is equal to 1 (the critical load) it means that in at least one state of the environment the total load over all classes must exceed 1, i.e. for at least one state we must have overload. This is true only for special cases of definition (i).

In this paper, special focus will be given to a multi-class single-server queue under the DPS discipline. The DPS discipline was first introduced by Kleinrock in [START_REF] Kleinrock | Time-shared systems: a theoretical treatment[END_REF] as an extension of the well-known egalitarian processor sharing (PS) discipline and has turned out to be very suitable to model the simultaneous parsing of diverse tasks, such as processing network data. Under this service discipline, the service capacity is divided between all present customers in proportion to their prescribed weights. Due to the challenging nature of DPS systems, most available results are in terms of limit theorems and moments. Fayolle et al. [START_REF] Fayolle | Sharing a processor among many job classes[END_REF] established the mean sojourn time conditioned on the service requirement. That analysis also yielded the mean queue lengths of the different classes, which were shown to depend on the entire service requirement distributions of all classes. The DPS model has finite mean queue lengths irrespective of any higher-order characteristics of the service distribution; see Avrachenkov et al. [START_REF] Avrachenkov | Discriminatory processor sharing revisited[END_REF]. This is an extension of a result for the Processor Sharing (PS) system, which holds while the queue is stable. DPS under a heavy-traffic regime was analysed in [START_REF] Grishechkin | On a relationship between processor sharing queues and Crump-Mode-Jagers branching processes[END_REF] assuming finite second moments of the service requirement distributions. Assuming exponential service requirement distributions, a direct approach to establish a heavytraffic limit for the joint queue-length distribution was described in [START_REF] Rege | Queue-length distribution for the discriminatory processor-sharing queue[END_REF] and extended to phase-type distributions in [START_REF] Verloop | Heavy-traffic analysis of a multiple-phase network with discriminatory processor sharing[END_REF]. Combining light-and heavy-traffic limits, in [START_REF] Izagirre | Interpolation approximations for the steady-state distribution in multi-class resource-sharing systems[END_REF] an interpolation approximation is derived for the steady-state distribution of the queue length and waiting time of DPS. The performance of DPS in overload is considered in [START_REF] Altman | DPS queues with stationary ergodic service times and the performance of TCP in overload[END_REF]. Asymptotics of the sojourn time have received attention in [START_REF] Borst | Tail asymptotics for discriminatory processor sharing with heavy-tailed service requirements[END_REF][START_REF] Borst | Sojourn time asymptotics in processor sharing queues[END_REF]. Game-theoretic aspects of DPS have been studied in [START_REF] Hassin | To Queue or Not to Queue: Equilibrium Behavior in Queueing Systems[END_REF][START_REF] Wu | Heavy traffic approximation of equilibria in resource sharing games[END_REF]. A thorough overview of DPS results can be found in [START_REF] Altman | A survey on discriminatory processor sharing[END_REF].

We are not aware of work analysing a DPS system under modulation. We refer to [START_REF] Núñez-Queija | Processor-sharing models for integrated-services networks[END_REF] where the processor sharing discipline (DPS discipline with unit weights) was analysed in a Markovian random environment. Multi-class queues in a random environment have been studied for different models in [START_REF] Budhiraja | Scheduling control for Markov-modulated single-server multiclass queueing systems in heavy traffic[END_REF][START_REF] Takine | Single-server queues with Markov-modulated arrivals and service speed[END_REF]. In [START_REF] Takine | Single-server queues with Markov-modulated arrivals and service speed[END_REF], a modulated system is studied where arrivals can only occur at transition epochs of the modulating process but service requirements are class dependent and generally distributed. Using a timechanging argument, the waiting time distribution is derived under the FCFS discipline. In [START_REF] Budhiraja | Scheduling control for Markov-modulated single-server multiclass queueing systems in heavy traffic[END_REF], the authors derive a Brownian control problem to establish a form of the cμ scheduling rule in heavy traffic under continuously modulated service requirements. By using a particular scaling, the time-scale separation of the external environment and the queue-length process is exploited. Similar scaling of a modulated queue can also be seen in results on the Markov-modulated infinite server queue in, for example, [START_REF] Blom | Markov-modulated infinite-server queues with general service times[END_REF].

The system we analyse in this paper is a single-server queue where the arrival rates, service requirements and service capacity are modulated. We focus on the setting where a customer's service requirement distribution is based on the state of the environment when it arrives and does not change throughout its service. The service capacity is, however, continuously modulated. This assumption is in line with the literature for various types of modulated queues; see [START_REF] Boxma | The M/G/1 queue with two service speeds[END_REF][START_REF] Dorsman | Heavy-traffic asymptotics for networks of parallel queues with Markov-modulated service speeds[END_REF][START_REF] Mahabhashyam | On queues with Markov modulated service rates[END_REF][START_REF] Takine | Single-server queues with Markov-modulated arrivals and service speed[END_REF]. In Remark 3 in Sect. 6, we discuss how part of our results can be extended to a more general model with continuously modulated service requirements. We derive the distribution of the workload under a heavy-traffic scaling for generally distributed service requirements and any service discipline which does not depend on the environment. We then turn our attention to the DPS discipline in a multi-class queue, which is a particular case of the general modulated system as described above. The weights of the DPS system, determining the service proportion, depend on a customer's class and do not change with the environment, which means that the workload result remains valid. An important finding in the present paper is that the queue-length vector under DPS becomes independent of the modulating process in the heavy-traffic limit, which is consistent with the modulated M/G/1 queue in, for example, [START_REF] Asmussen | The heavy traffic limit of a class of Markovian queueing models[END_REF][START_REF] Dimitrov | Single-server queueing system with Markov-modulated arrivals and service times[END_REF]. This, together with the obtained result on the workload, allows us to derive the distribution of the queue-length vector under DPS and to show that it undergoes a state-space collapse.

The remainder of the paper is organized as follows. In Sect. 2 we describe the model. In Sect. 3 we study the workload of a single-server queue with modulated arrivals, service capacity and service requirement distribution, establishing its distribution in heavy traffic. From Sect. 4 onwards the focus is on DPS. In Sect. 4 we derive some basic properties of the queue-length distribution, obtain a rate conservation law and derive an equation for the moments of the queue lengths weighted with the modulated service capacity. Section 5 is devoted to the heavy-traffic scaling; there we show that the distribution of the environment becomes independent of that of the queuelength vector, in addition to deriving two technical lemmas. The exponential limiting distribution of the joint queue length in heavy traffic follows in Sect. 6. The result is shown in two steps; Sect. 6.1 for the state-space collapse and Sect. 6.2 for the exact limiting distribution, where we rely on the workload result of Sect. 3. We conclude with a summary and some open questions in Sect. 7.

Model

We analyse a single-server queue modulated by an independent external environment, which is formalized by an irreducible continuous time Markov chain on a finite statespace {1, . . . , D}. The modulating process is denoted by Z and is governed by an infinitesimal generator matrix Q = (q d ) D d, =1 with an invariant distribution π = (π 1 , . . . , π D ). In what follows, vectors are generally denoted in bold. New customers arrive according to a Poisson distribution with rate λ d when the environment is in state d. A customer arriving in state d has a service requirement distribution given by a function H d (•) with Laplace-Stieltjes transform (LST) h d (•). The service requirement does not change further with the environment. The first and second moments are given by h d1 and h d2 , respectively. In addition we let the service capacity be scaled by a factor c d during the environment's stay in state d; this can thus change during the service of the customers. The traffic intensity will be measured as

ρ ∞ = d π d λ d h d1 /c ∞ , (1)
with c ∞ := d π d c d being the service capacity averaged over the environment. The workload is defined as the time it takes to empty the system at an arbitrary moment in time given the observed environment and is denoted by W . In Sect. 3 we study the workload and the environment as a two-dimensional process (W, Z), under any service discipline that is independent of the environment.

The first main result of this paper concerns the distribution of the workload when the traffic intensity approaches its critical point. The system is said to be in heavy traffic when ρ ∞ approaches 1. Let N > 0 and define the following parametrization:

λ (N ) d := λ d ρ ∞ (1 -1/N ) → λ d ρ ∞ =: λd , as N → ∞, (2)
where ρ ∞ is based on the unscaled parameters. Prelimit quantities will be denoted with a superscript (N ); the prelimit traffic intensity is thus ρ

(N )
∞ and is equal to 1 -1/N . Limiting quantities will have a ˆ; in heavy traffic the traffic intensity is denoted ρ∞ and is equal to 1.

In the remainder of the paper, starting in Sect. 4, we analyse a single-server queue with K customer classes under the discriminatory processor sharing policy; the queue is again embedded in a random environment. Let α k,d be the probability that a customer, arriving while the environment is in state d, is of class k; note that k α k,d = 1 for a given d. The Poisson arrival rate of a class-k customer is denoted by λ k,d := α k,d λ d and for each class k it is assumed that λ k,d > 0 for at least one state d. In the multiclass setting we assume that a class-k customer has an exponentially distributed service requirement with mean 1/μ k . We believe that the results obtained in this paper can be extended to phase-type distributed service requirements, the latter being dense in the space of all distributions on [0, ∞). For the non-modulated DPS queue, the phase-type analysis was performed in [START_REF] Verloop | Heavy-traffic analysis of a multiple-phase network with discriminatory processor sharing[END_REF] using similar proof techniques. For ease of exposition, however, we focus here on the exponential case.

We no longer let the service requirement of a particular customer be environment dependent (although the distribution of an arbitrary customer is, as explained in Sect. 6.2). One can, however, retrieve the environment-dependent service requirements by introducing additional classes for each environment; see Remark 1. By referring to a class-k customer's service rate while in state d as μ k,d := μ k c d , we take the modulated service capacity into account. Most of the results for the queue length can in fact be shown without assuming this product form, representing a system where the service requirements are continuously modulated; see Remark 3, Sect. 6.1, for further details.

We denote the average arrival rate of class-k customers by λ k,∞ := d λ k,d π d ; similarly we denote the average service rate for class

k by μ k,∞ := d μ k,d π d = μ k c ∞ and ρ k,∞ := λ k,∞ /μ k,∞ .
The aggregate traffic intensity for the multi-class model is defined as

ρ ∞ := K k=1 ρ k,∞ ,
which is consistent with the definition in Eq. [START_REF] Altman | DPS queues with stationary ergodic service times and the performance of TCP in overload[END_REF].

Let the state of the multi-class system be described by the vector of random variables (M 1 , . . . , M K , Z ) =: (M, Z ), where M k is the number of class-k customers, for k = 1, . . . , K . As before, Z represents the state of the background process. In a DPS system, the random fraction of service given to a class-k customer is

d = 1 d = 2 d ∈ {1, 2} (1,1) λ 1,1 0 μ 1,1 g 1 (1,2) 0 λ 1,2 μ 1,2 g 1 (2,1) λ 2,1 0 μ 2,1 g 2 (2,2) 0 λ 2,2 μ 2,2 g 2 g k j g j M j
, where g k are weight parameters associated with each class k.

When in heavy traffic, we denote ρk,∞ := λk,∞ /μ k,∞ and thus have for the multi- From Table 1 one can easily see how a K = D = 2 system can be written into one with K = 4 and D = 2. The arrival rates are still modulated, but in an on-off way. The service rates μ k,d are now non-modulated.

class system k ρk,∞ = k λk,∞ μ k,∞ = 1 ρ ∞ k λ k,∞ μ k,∞ = 1 ρ ∞ k ρ k,∞ = 1.

Workload

In this section we consider the workload in a modulated queue and extend the results for an M/G/1-type queue from Falin and Falin [START_REF] Falin | Heavy traffic analysis of M/G/1 type queueing systems with Markov-modulated arrivals[END_REF] and Dimitrov [START_REF] Dimitrov | Single-server queueing system with Markov-modulated arrivals and service times[END_REF] to include modulated service capacity. We derive the mean of the workload and then its distribution in the heavy-traffic regime.

Let p 0,d = P(W = 0, Z = d) and let a = (a 1 , . . . , a D ) T be a vector solving

[Q • a] d = c d -λ d h d1 -c ∞ (1 -ρ ∞ ), (3)
for d = 1,..., D. Note that such a solution always exists since the right-hand side vector of Eq. ( 3) is orthogonal to π. We obtain the following result for the mean workload.

Proposition 3.1 For any service requirement distribution H d (•) and any service discipline that does not depend on the state of the environment, the mean of the workload satisfies

EW = d π d λ d h d2 /2 + a d π d (λ d h d1 -c d ) + p 0,d c d a d c ∞ (1 -ρ ∞ ) , ( 4
)
where a is a solution of Eq. (3). Furthermore,

1 -ρ ∞ = d p 0,d c d /c ∞ .
Remark 2 Although the solution vector a is not unique, the term

D d=1 a d [π d (λ d h d1 -c d ) + p 0,d c d ],
as appearing in Eq. ( 4), is. This is due to the following argument: Suppose a and a * are two solutions to Eq. ( 3). Then 0 = Qa -

Qa * = Q(a -a * ), so (a -a * ) is in the nullspace of Q.
But Q is a generator matrix so it can easily be seen that Q r = 0 for any vector r • 1 T = (r, r, . . . , r ) T , r ∈ R. Also, since the environment is an irreducible Markov chain, the nullspace of Q has dimension 1, and therefore, (a

-a * ) = r a • 1 T , for some r a ∈ R. Thus, D d=1 (a d -a * d )[π d (λ d h d1 -c d ) + p 0,d c d ] = r a c ∞ (ρ ∞ -1) + r a c ∞ (1 -ρ ∞ ) = 0,
where the first term follows by the definition of ρ ∞ and c ∞ and the second term comes from Eq. ( 7) in the proof below.

Proof of Proposition 3.

F d (x, t + dt) = (1 -λ d dt + q dd dt)F d (x + c d dt, t) + =d q d F (x + c dt, t)dt + λ d dt x 0 F d (x + c d dt -y, t)dH d (y).
We let t → ∞ to go to steady state, and since

F d (x + c d dt) -F d (x) dt = c d F d (x + c d dt) -F d (x) c d dt -→ dt↓0 c d F d (x),
we obtain

c d F d (x) = (λ d -q dd )F d (x) - =d q d F (x) -λ d x 0 F d (x -y)dH d (y). Denote the LST of F d (•) by ϕ d (s) = E[e -sW (t) , Z (t) = d] = p 0,d + ∞ 0+ e -sx dF d (x). The corresponding transform equation becomes D =1 q d ϕ (s) = [λ d (1 -h d (s)) -sc d ]ϕ d (s) + sp 0,d c d .
(

It is now convenient to sum over d and divide through Eq. ( 5) with s to get zero on the left-hand side, leading to

d p 0,d c d = d ϕ d (s) c d -λ d (1 -h d (s)) s . ( 6
)
Using that ϕ d (0) = π d and by l'Hôpital's rule

lim s↓0 (1 -h d (s)) s = -lim s↓0 h d (s) = h d1 .
We get by taking the limit s → 0 of Eq. ( 6) that

d p 0,d c d c ∞ = 1 -ρ ∞ . (7)
We differentiate Eq. ( 6) w.r.t. s:

d ϕ d (s)λ d h d (s) s + 1 -h d (s) s 2 = d ϕ d (s) c d -λ d 1 -h d (s) s , which in the limit s → 0 results in d π d λ d h d2 2 = d W d [c d -λ d h d1 ] , (8)
with the first moment of the workload while in state d being

-lim s↓0 ϕ d (s) = E[W, Z = d] =: W d .
Now multiply Eq. ( 5) with a d , sum over d, take the derivative w.r.t. s and let s → 0 to obtain

- d W d [c ∞ (ρ ∞ -1) + c d -λ d h d1 ] = d a d π d (λ d h d1 -c d ) + p 0,d c d a d ,
by using Eq. [START_REF] Asmussen | The heavy traffic limit of a class of Markovian queueing models[END_REF]. Adding this equation to Eq. ( 8) yields

c ∞ (1 -ρ ∞ ) d W d = d π d λ d h d2 2 + a d π d (λ d h d1 -c d ) + p 0,d c d a d , (9)
which gives the desired expression for the mean of the workload, EW = d W d .

Eq. ( 7) makes it clear that in heavy traffic, that is, when ρ (N )

∞ → 1, all the probabilities p (N ) 0,d go to zero. Also, in heavy traffic, the right-hand side of Eq. ( 3) reduces to c d -λd h d1 . Recalling the parametrization 1 -ρ (N ) ∞ = 1/N in Sect. 2, we obtain the following result.

Proposition 3.2 For any service distribution H d (•) and any service discipline that does not depend on the state of the environment, the mean of the workload in heavy traffic satisfies

lim N →∞ 1 N EW (N ) = 1 c ∞ d π d λd h d2 /2 + a d ( λd h d1 -c d ) , ( 10
)
where a is a solution of

[Q • a] d = c d -λd h d1 .
Proof Under the heavy-traffic scaling, the empty probabilities p (N ) 0,d go to zero for d = 1, . . . , D. The result then follows immediately from Eq. ( 4) and 1-ρ

(N ) ∞ = 1/N .
This leads to the main result of this section.

Theorem 3.3

In heavy traffic, the scaled workload 1 N W (N ) converges in distribution to Ŵ , where Ŵ is exponentially distributed with mean given in Eq. [START_REF] Boxma | The M/G/1 queue with two service speeds[END_REF].

Proof This follows from combining Proposition 3.2 with Theorem 4 in [START_REF] Dimitrov | Single-server queueing system with Markov-modulated arrivals and service times[END_REF]. The full proof is in the Appendix.

The above results yields that Ŵ is relatively compact, which together with the metric space being separable and complete implies that the scaled workload 1 N W (N ) is tight, by Prohorov's theorem [START_REF] Billingsley | Convergence of Probability Measures[END_REF].

Queue-length vector under DPS

In the remainder of the paper we focus on the multi-class model under DPS, where we assume that the service requirements of class-k customers are exponentially distributed with rate μ k . In this section we establish some properties of the joint queue-length distribution. We start with the flow equations, followed by a rate conservation law and an equation for the moments of the queue lengths conditioned on the environment.

Equating the flow in and out of state

(M, Z ) = (m, d) yields (noting that -q dd = =d q d ) K k=1 λ k,d + g k m k i g i m i μ k,d • 1 {m k >0} -q dd p m,d = K k=1 λ k,d p m-e k ,d • 1 {m k >0} + g k (m k + 1) i g i m i + g k μ k,d p m+e k ,d + =d q d p m, , (11)
where p m,d := P((M, Z ) = (m, d)) and e k is the vector with 1 in the k-th place and zeros elsewhere. We now define the partial probability generating function (PGF) for when the background process is in state d:

P d (z) := E[z M 1 1 • • • z M K K • 1 {Z =d} ] := ∞ m 1 =0 • • • ∞ m K =0 P(M 1 = m 1 , . . . , M K = m K , Z = d) • z m 1 1 • • • z m K K = m≥0 p m,
) := E[z M 1 1 • • • z M K K ] = D d=1 P d (z).
We also define

R d (z) := m≥0 p m,d z m j g j m j • 1 K j=1 m j >0 , hence ∂ R d (z) ∂z k = z -1 k m≥e k m k j g j m j p m,d z m • 1 K j=1 m j >0
.

By multiplying the flow equation [START_REF] Budhiraja | Scheduling control for Markov-modulated single-server multiclass queueing systems in heavy traffic[END_REF] with z m , summing over all vectors m ≥ 0 and rearranging terms, we can eventually write it in terms of the PGF P d (z) and the partial derivative

∂ R d (z)/∂z k , that is, K k=1 λ k,d (1 -z k )P d (z) + μ k,d g k (z k -1) ∂ R d (z) ∂z k = D =1 P (z)q d . ( 12
)
It will be convenient to write the equation fully in terms of ∂ R d /∂z k , so we note the relation

P d (z) = K k=1 g k z k ∂ R d (z) ∂z k + p 0,d , (13)
where p 0,d = P((M, Z ) = (0, d)) is the probability of an empty queue in state d, which is equivalent to the probability of no workload defined in Sect. 3. We incorporate this into Eq. ( 12) to obtain

K k=1 λ k,d (1 -z k ) ⎡ ⎣ K j=1 g j z j ∂ R d (z) ∂z j + p 0,d ⎤ ⎦ + K k=1 μ k,d g k (z k -1) ∂ R d (z) ∂z k = D =1 K k=1 g k z k ∂ R (z) ∂z k + p 0, q d . ( 14
)
This equation we will use later when deriving the heavy-traffic limit.

For the M/M/1 queue with modulated arrivals and service times, moments of the queue length and the sojourn times can be found for the FCFS service discipline in [START_REF] Yechiali | Queueing problems with heterogeneous arrival and service[END_REF] and [START_REF] Mahabhashyam | On queues with Markov modulated service rates[END_REF], respectively. In [START_REF] Rege | Queue-length distribution for the discriminatory processor-sharing queue[END_REF] the authors establish a recursive formula to calculate moments of the queue length in a non-modulated DPS system. In a similar fashion, we obtain an expression for the sum of the state-dependent moments weighted with the capacity of the server. We also derive a rate conservation law, which shows how the average arrival rates per class are proportional to the resources allocated to that same class, and the service they receive. Both results can be found in the following proposition.

Proposition 4.1 When the queue is stable, the average number of class-k arrivals is proportional to the service resources allocated to class-k customers, i.e.

λ k,∞ = d μ k,d E g k M k j g j M j • 1 { j M j >0} • 1 {Z =d} , for k = 1, . . . , K . Furthermore, the state-dependent expectations of M k satisfy d c d E[M k • 1 {Z =d} ] = λ k,∞ μ k + d, j g j λ k,d E[M j • 1 {Z =d} ] + λ j,d E[M k • 1 {Z =d} ] μ k g k + μ j g j .
Proof We sum Eq. ( 12) over d and take the derivative w.r.t. z i to obtain

d ⎡ ⎣ -λ i,d P d (z) + j λ j,d (1 -z j ) ∂ P d (z) ∂z j + μ i,d g i ∂ R d (z) ∂z i + j μ j,d g j (z j -1) ∂ 2 R d (z) ∂z i ∂z j ⎤ ⎦ = 0.
Letting z → 1 yields

d μ i,d g i ∂ R d (z) ∂z i z→1 -λ i,d P d (z) z→1 = 0. ( 15
)
Since

lim z→1 P d (z) = π d ,
the following conservation law results from Eq. ( 15), for k = 1, . . . , K :

λ k,∞ = d μ k,d g k ∂ R d (z) ∂z k z→1 = d μ k,d E g k M k j g j M j • 1 { j M j >0} • 1 {Z =d} . ( 16
)
By taking partial derivatives of Eq. ( 13), we obtain after standard calculations the recursive relation

∂ j P d (z) ∂z i 1 • • • ∂z i j z→1 = K k=1 g k ∂ j+1 R d (z) ∂z i 1 • • • ∂z i j ∂z k + j =1 g i ∂ j R d (z) ∂z i 1 • • • ∂z i j . ( 17
)
In particular, this yields the explicit form

E[M k • 1 {Z =d} ] = ∂ P d ∂z k z→1 = j g j ∂ 2 R d ∂z k ∂z j z→1 + g k ∂ R d ∂z k z→1 .
Proceeding from the rate conservation law, Eq. ( 16), and by using μ k,d = μ k c d , we have

d c d g k ∂ R d ∂z k z→1 = λ k,∞ μ k .
By taking two partial derivatives of the balance equation Eq. ( 12), we can solve for a second mixed derivative of R d , namely

∂ 2 R d ∂z k ∂z j z→1 = E M k M j •1 {Z = } q d +λ k,d E M j •1 {Z =d} +λ j,d E M k •1 {Z =d} μ k,d g k +μ j,d g j ,
thus also yielding a mixed moment. Summing over the weighted moments of the number of class-k customers while in state d, we obtain a linear equation resembling Eq. ( 16) in [START_REF] Rege | Queue-length distribution for the discriminatory processor-sharing queue[END_REF] for the non-modulated DPS queue:

d c d E[M k • 1 {Z =d} ] = λ k,∞ μ k + d, j c d g j ∂ 2 R d ∂z k ∂z j z→1
(by Eq. ( 17))

= λ k,∞ μ k + d, j g j E[M k M j •1 {Z = } ]q d +λ k,d E[M j •1 {Z =d} ]+λ j,d E[M k •1 {Z =d} ] μ k g k +μ j g j = λ k,∞ μ k + d, j g j λ k,d E[M j • 1 {Z =d} ] + λ j,d E[M k • 1 {Z =d} ] μ k g k + μ j g j ,
where the last equality comes from d q d = 0.

Preliminary results for the queue length in heavy traffic

We proceed to show that in heavy traffic the distribution of the environment and the joint queue length become independent. This result, along with two technical lemmas that we derive in this section, will later help establish the main result about the limiting queue length under DPS, presented in Sect. 6. Here we consider the queue-length vector (M 1 , . . . , M K ) scaled with 1/N and evaluate the PGF in z 1/N . The objective is to determine the distribution of 1

N (M (N ) 1 , . . . , M (N ) K )
• 1 {Z =d} as N goes to infinity. We will state the existence of the limiting vector, and thus also the limit of the generating function P (N ) d (z 1/N ), as an assumption. This assumption will be proven in Sect. 6.2. The superscript N denotes dependency on the prelimit parameter λ (N ) d . We make use of the change of variables e -s k = z k , for s k > 0, and denote e -s/N := (e -s 1 /N , . . . , e -s K /N ). Assuming that the limit exists, we use the new variables to define the heavy-traffic quantities: We let p0,d := lim N →∞ p (N ) 0,d , Pd (s) := lim N →∞ P (N ) d (e -s/N ) = lim N →∞ E[e -j s j M j /N • 1 {Z =d} ] and P(s) := d Pd (s) = lim N →∞ E[e -j s j M j /N ]. We denote by ( M1 , . . . , MK ) the random vector corresponding to the LST P(s). Finally, let

Rd (s) := E 1 -e -j s j M j j g j M j • 1 { j M j >0} • 1 {Z =d} ,
where the 1 in the numerator is to ensure that the bracketed expression remains bounded when the queue length quantities Mk are all near zero. We can now proceed to the following lemma.

Lemma 5.1 If lim N →∞ P (N )

d (e -s/N ) exists, then it satisfies

Pd (s) = K k=1 g k ∂ Rd (s) ∂s k + p0,d . ( 18
)
Proof From Eq. ( 13), lim

N →∞ P (N ) d (e -s/N ) = lim N →∞ K k=1 g k z k ∂ R (N ) d (z) ∂z k z=e -s/N + p0,d . ( 19
)
Note that lim

N →∞ z k ∂ R (N ) d (z) ∂z z=e -s/N = lim N →∞ m≥e k m k j g j m j p (N ) m,d z m • 1 { j m j >0} z=e -s/N = lim N →∞ m≥e k m k j g j m j p (N ) m,d e -s 1 m 1 /N • • • • • e -s K m K /N • 1 { j m j >0} = lim N →∞ E M (N ) k j g j M (N ) j e -j s j M (N ) j /N • 1 { j M (N ) j /N >0} • 1 {Z =d} = E Mk j g j M j e -j s j M j • 1 { j M j >0} • 1 {Z =d} = ∂ Rd (s) ∂s k . ( 20
)
The second-to-last step follows from the fact that

M (N ) k j g j M (N ) j • e -j s j M (N ) j • 1 { j M (N ) j >0}
is upper bounded by 1/ min j (g j ). By the continuous mapping theorem (see Billingsley [START_REF] Billingsley | Convergence of Probability Measures[END_REF]), it converges in distribution to

Mk j g j M j • e -j s j M j • 1 { j M j >0} .
The environment is not affected by the heavy-traffic scaling. Eqs. ( 19) and ( 20) now conclude the proof.

With the help of Lemma 5.1 we obtain: Proposition 5.2 If lim N →∞ P (N ) d (e -s/N ) exists, the joint queue-length distribution is independent of the environment in heavy traffic, that is,

Pd (s) = π d P(s).
Proof We use the change of variables z k = e -s k . Since

z k ∂ R d ∂z k (z) z=e -s = - ∂ R d ∂s k (e -s ),
we obtain, by applying the heavy-traffic scaling to Eq. ( 14),

k λ (N ) k,d 1 -e -s k /N ⎡ ⎣ K j=1 g j -∂ R (N ) d ∂s j e -s/N + p (N ) 0,d ⎤ ⎦ -μ k,d g k e -s k /N -1 e s k /N ∂ R (N ) d ∂s k e -s/N = D =1 K k=1 g k -∂ R (N ) ∂s k e + p (N ) 0, q d .
With Taylor expansion we obtain

k λ (N ) k,d s k N - s 2 k N 2 ⎡ ⎣ K j=1 g j -∂ R (N ) d ∂s j e -s/N + p (N ) 0,d ⎤ ⎦ -μ k,d g k s k N + s 2 k N 2 -∂ R (N ) d ∂s k e -s/N = D =1 K k=1 g k -∂ R (N ) ∂s k e -s/N + p (N ) 0, q d + O(N -3 ). ( 21
)
Since

-∂ R (N ) d ∂s j
is bounded (see the proof of Lemma 5.1) and converges to ∂ Rd ∂s j , we obtain, as

N → ∞, ν • [Q] d = D =1 K k=1 g k ∂ R (s) ∂s k + p0, q d = 0, s ≥ 0, ∀d, (22)
where [Q] d is the d th column of Q and ν is a row vector with ν = K k=1 g k ∂ R (s) ∂s k + p0, . This implies that ν Q = 0, and since Q is a generator we conclude that

ν d = K k=1 g k ∂ Rd (s) ∂s k + p0,d = π d x,
where x does not depend on d. Observe now that by Lemma 5.1 we have

Pd (s) -p0,d = K k=1 g k ∂ Rd (s) ∂s k = π d x -p0,d = E[1 {Z =d} ]x -p0,d .
Since Pd (s) = E e -j s j M j • 1 {Z =d} , this implies that x = E e -j s j M j = P(s). This shows that the environment becomes independent from the joint queuelength process in the heavy-traffic limit.

The flow equation, Eq. ( 14), simplifies considerably in heavy traffic, as shown in the following lemma.

Lemma 5.3 If lim N →∞ P (N )

d (e -s/N ) exists, then Rd (s) satisfies the following equation:

0 = k,d F k,d (s) ∂ Rd (s) ∂s k , ∀s ≥ 0, with F k,d (s) defined as F k,d (s) := g k ⎛ ⎝ j λ j,d s j -μ k,d s k ⎞ ⎠ .
Proof We start by multiplying through Eq. ( 21) with N , followed by summing over d. Due to Q being a generator, this eliminates the right-hand side with the transition rates

q d : k,d ⎡ ⎣ λ (N ) k,d s k - s 2 k N ⎡ ⎣ K j=1 g j -∂ R (N ) d ∂s j (e -s/N ) + p (N ) 0,d ⎤ ⎦ ⎤ ⎦ - k,d μ k,d g k s k + s 2 k N -∂ R (N ) d ∂s k (e -s/N ) + O(N -2 ) = 0. Taking the limit N → ∞ yields 0 = k,d λk,d s k j g j ∂ Rd (s) ∂s j - k,d μ k,d g k s k ∂ Rd (s) ∂s k = k,d g k ∂ Rd (s) ∂s k j λ j,d s j - k,d μ k,d g k s k ∂ Rd (s) ∂s k = k,d g k ⎛ ⎝ j λ j,d s j -μ k,d s k ⎞ ⎠ ∂ Rd (s) ∂s k = k,d F k,d (s) ∂ Rd (s) ∂s k . ( 23
)
In what follows we focus on

F k,∞ (s) := d F k,d (s)π d = g k ⎛ ⎝ j λ j,∞ s j -μ k,∞ s k ⎞ ⎠ ,
and denote its vector counterpart by F ∞ (s) = (F 1,∞ (s), . . . , F K ,∞ (s)).

Queue-length distribution in heavy traffic

We now state and subsequently prove our main result about the queue-length distribution.

Theorem 6.1 When scaled by

1/N = (1 -ρ (N )
∞ ), the queue-length vector converges in distribution as (ρ

(N ) 1,∞ , . . . , ρ (N ) K ,∞ ) → ( ρ1,∞ , . . . , ρK,∞ ), i.e. ρ (N ) ∞ → 1, namely 1 N M (N ) 1 , . . . , M (N ) K • 1 {Z =d} d → M1 , . . . , MK • 1 {Z =d} d = π d ρ1,∞ g 1 , . . . , ρK,∞ g K • X, ( 24
)
where d → denotes convergence in distribution and X is exponentially distributed with mean

EX = k ρk,∞ /μ k -d c d π d a d (1 -ρd ) c ∞ k ρk,∞ /(g k μ k ) , ( 25
)
with ρd := c -1 d k λk,d /μ k and a = (a 1 , . . . , a D ) T being a solution of [Q • a] d = c d (1 -ρd ).
We will prove this theorem in the two following subsections, first showing in Sect. 6.1 the state-space collapse observed in Eq. ( 24) and then in Sect. 6.2 we will show that X is exponentially distributed with the mean given by Eq. (25).

State-space collapse

The first part of the proof of Theorem 6.1 is the state-space collapse. In this section, we assume lim

N →∞ P (N )
d (e -s/N ) exists. Observe that due to Proposition 5.2,

Rd (s) = E 1 -e -k s k Mk k g k Mk • 1 { k Mk >0} • 1 {Z =d} = E 1 -e -k s k Mk k g k Mk • 1 { k Mk >0} • π d = : R(s)π d , (26)
where the last equation defines R(s), which is independent of d. We now derive some properties of R(s). Lemma 6.2 R(s) is constant on a (K -1)-dimensional hyperplane H c , where

H c := s ≥ 0 : k ρk,∞ g k s k = c , c > 0.
Proof We follow closely the steps of the proof of Lemma three in [START_REF] Verloop | Heavy-traffic analysis of a multiple-phase network with discriminatory processor sharing[END_REF]. The proof has three steps: (i) Show that F k,∞ (s) is parallel to the hyperplane. Hence, any flow corresponding to F k,∞ that starts in the plane, stays in the plane. (ii) Show that R(s) is constant along each flow in the hyperplane and (iii) show that each flow in the hyperplane converges to a unique point. This implies that R(s) is constant on the hyperplane.

(i) F ∞ (s) is parallel to H c Observe that with 1 = k ρk,∞ and ρk,∞ = λk,∞ /μ k,∞ , k ρk,∞ g k F k,∞ (s) = k ρk,∞ ⎛ ⎝ j λ j,∞ s j -μ k,∞ s k ⎞ ⎠ = j λ j,∞ s j - k ρk,∞ μ k,∞ s k = k λk,∞ s k - k λk,∞ s k = 0.
This indicates that the K -dimensional vector F ∞ (s) is parallel to the hyperplane. (ii) R(s) is constant along flows in H c For each state s ≥ 0, there exists a unique flow

f (u) = ( f 1 (u), . . . , f K (u)) T parametrized by u ≥ 0, such that f (0) = s and d f k (u) du = F k,∞ ( f (u)). ( 27
)
Due to (i), any flow that starts in H c , stays in

H c . Now, d R( f (u)) du = K k=1 d f k (u) du • ∂ R(s) ∂s k s= f (u) = K k=1 F k,∞ ( f (u)) • ∂ R(s) ∂s k s= f (u) = K k=1 D d=1 F k,d ( f (u))π d • ∂ R(s) ∂s k s= f (u) = K k=1 D d=1 F k,d ( f (u)) • ∂ Rd (s) ∂s k s= f (u)
= 0, by Eq. ( 23),

implying that R( f (u)) is constant along each flow f (u) which lies in H c . (iii) Each flow in H c converges to a unique point
Here we first write the flow specifications in a vector-matrix form, then show that one eigenvalue of that matrix is zero with eigenvector s * ∈ H 1 , and the other eigenvalues are negative, and thus we can write

f (u) = c • s * + g(u)
where lim u→∞ g(u) = 0. Eq. ( 27) can be written in matrix-vector form as

f (u) = A f (u), with A = ⎛ ⎜ ⎜ ⎜ ⎝ g 1 ( λ1,∞ -μ 1,∞ ) g 1 λ2,∞ • • • g 1 λK,∞ g 2 λ1,∞ g 2 ( λ2,∞ -μ 2,∞ ) • • • g 2 λK,∞ . . . . . . . . . g K λ1,∞ • • • g K ( λK,∞ -μ K ,∞ ) ⎞ ⎟ ⎟ ⎟ ⎠ .
Let D be the diagonal matrix with d i = ρi,∞ /g i on the diagonal. Then with

S:=D AD -1 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ g 1 ( λ1,∞ -μ 1,∞ ) g 2 ρ1,∞ ρ2,∞ λ2,∞ • • • g K ρ1,∞ ρK,∞ λK,∞ g 1 ρ2,∞ ρ1,∞ λ1,∞ g 2 ( λ2,∞ -μ 2,∞ ) • • • g K ρ2,∞ ρK,∞ λK,∞ . . . . . . . . . g 1 ρK,∞ ρ1,∞ λ1,∞ • • • g K ( λK,∞ -μ K ,∞ ) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .
S T is a generator corresponding to a finite-state Markov chain.

From the proof of Lemma 4 in [START_REF] Verloop | Heavy-traffic analysis of a multiple-phase network with discriminatory processor sharing[END_REF], it is easily seen that the Markov chain corresponding to S T is irreducible (since we assume that all λ k,d > 0 for all k and at least one d). Retracing the arguments stated there, for completeness, it follows that this Markov chain has a unique equilibrium distribution (column) vector, η, such that η T S T = 0. In particular, 0 is an eigenvalue with multiplicity one and all other eigenvalues have a strictly negative real part; see [START_REF] Asmussen | Applied Probability and Queues[END_REF]. Since the eigenvalues of S T and A are the same, 0 is also an eigenvalue of A with corresponding right eigenvector s * = D -1 η, s * ≥ 0, s * ∈ H 1 . The solution of the linear system f (u) = A f (u), f (0) ∈ H c can now be written as the sum of the homogeneous and the particular solution, i.e. f (u) = c • s * + g(u), where lim u→∞ g(u) = 0. This implies that all the flows in H c converge to one common point c • s * .

Combining (i), (ii) and (iii), we conclude that the function R(s) is constant on H c .

As a consequence of Lemma 6.2, the function R(s) depends on s only through the sum K k=1 ( ρk,∞ /g k )s k . Therefore, there exists a function R * : R → R such that

R(s) = R * ( K k=1 ( ρk,∞ /g k )s k ). Then ∂ ∂s k R(s) = ρk,∞ g k d R * (v) dv v= K k=1 ( ρk,∞ /g k )s k , so we obtain E[e -K k=1 s k Mk ] = lim N →∞ D d=1 P (N ) d (e -s/N ) = D d=1 K k=1 g k ∂ Rd (s) ∂s k = D d=1 K k=1 g k ∂ R(s) ∂s k π d = K k=1 ρk,∞ d R * (v) dv v= K k=1 ( ρk,∞ /g k )s k = d R * (v) dv v= K k=1 ( ρk,∞ /g k )s k , ( 28
)
which only depends on v = K k=1 ( ρk,∞ /g k )s k . Since we also have

E[e -K k=1 s k Mk ] = E ⎡ ⎣ e - g 1 ρ1,∞ v M1 • e -ρ2,∞ g 2 s 2 g 2 ρ2,∞ M2 - g 1 ρ1,∞ M1 • • • • • e -ρK ,∞ g K s K g K ρK ,∞ MK - g 1 ρ1,∞ M1 ⎤ ⎦ ,
this, together with Eq. ( 28), implies that

g j ρ j,∞ M j -g 1 ρ1,∞ M1 = 0, for all j = 1,..., K . Thus (g k /ρ ˆk,∞ )M ˆ k d = (g j /ρ ˆ j,∞ )M ˆ j ,
for all k, j, almost surely. Combining this finding with that of Eq. ( 22), we obtain Eq. ( 24) with X distributed as

(g 1 /ρ ˆ1,∞ )M ˆ 1 .
Remark 3 (Continuously modulated service requirements) In Sect. 3 we saw that the critical load is indeed reached when ρ ∞ → 1, since then p 0,d → 0. This indicates that (1 -ρ ∞ ) is the right heavy-traffic scaling when μ k,d = μ k c d .Thisisless clear for a general μ k,d , that is, for continuously modulated service requirements, where the environment can influence the departure rate of customers present in the system. For that setting, the workload process is no longer independent of the employed scheduling discipline, since the decision on which class to serve impacts the rate at which customers leave. We are not aware of any results on workload and waiting time distributions where the service distribution is a general function of both class and environment.

The majority of the preceding queue-length results in this paper can, however, be proven without the restriction of the product form, i.e. for continuously modulated service requirements. The traffic intensity for this variant is defined as for the multiclass model above, only this time one cannot split the average class-k service rate into μ k,∞ = μ k c ∞ . The traffic intensity per class k, ρ k,∞ = λ k,∞ /μ k,∞ , is in line with the Markov-modulated single-server queues. Assuming there exists a scaling f (N ) such that f (N )(M 1 , . . . , M K )1 {Z =d} converges in distribution, it can be shown that the empty probabilities p 0,d vanish in heavy traffic as N → ∞, for a general μ k,d . This property then follows from Proposition 5.2 without relying on the workload results from Sect. 3 and the product form assumed there. Furthermore, under this assumption, all results in Sect. 6.1 hold, implying that a state-space collapse will appear. In other words, we can prove the first half of Theorem 6.1. However, we do not know what the distribution of the common factor X will be.

Distribution of the common factor

In order to prove that the limiting queue-length distribution exists and to find the common factor of the queue-length distribution in heavy traffic, the random variable X , we make use of the results on the workload of the total system. From [START_REF] Verloop | Heavy-traffic analysis of a multiple-phase network with discriminatory processor sharing[END_REF] and Eq. [START_REF] Kleinrock | Time-shared systems: a theoretical treatment[END_REF] we have that

Ŵ d = K k=1 Mk μ k = X • K k=1 ρk,∞ g k μ k . ( 29
)
In order to apply the workload result of Sect.

H d (x) := K k=1 α k,d H k (x).
The overall service requirement distribution thus depends on the state of the environment at its arrival. With exponential service requirements, the corresponding LST is given by

h d (s) = K k=1 α k,d μ k μ k + s , s ≥ 0, (30)
and the first and second moment are given by

h d1 = k α k,d μ k , h d2 = 2 k α k,d μ 2 k . ( 31
)
We can now apply the result of 3.3 with moments as given in Eq. [START_REF] Takine | Single-server queues with Markov-modulated arrivals and service speed[END_REF]. Hence, we have that Ŵ is exponentially distributed with mean

E Ŵ = c -1 ∞ k ρk,∞ /μ k - d c d π d a d (1 -ρd ) ,
where a is a solution of

[Q • a] d = c d -λd k α k,d μ k = c d (1 -ρd ).
Along with Eq. ( 29), this yields the mean of the exponential random variable X :

EX = E Ŵ k ρk,∞ /(g k μ k ) = k ρk,∞ /μ k -d c d π d a d (1 -ρd ) c ∞ k ρk,∞ /(g k μ k ) . ( 32
)
The first term of the numerator is in accordance with the results of [START_REF] Rege | Queue-length distribution for the discriminatory processor-sharing queue[END_REF] and [START_REF] Verloop | Heavy-traffic analysis of a multiple-phase network with discriminatory processor sharing[END_REF], and the second term is a result of the random environment.

The results in Sects. 5 and 6.1 are based on the assumption that lim N →∞ 1 N M • 1 {Z =d} exists. Since the scaled workload is tight, see Sect. 3, so is the scaled queue length. Then, by Prohorov's theorem ( [START_REF] Billingsley | Convergence of Probability Measures[END_REF]) there exists a subsequence of N such that 1 N M k converges in distribution, and hence for this subsequence lim N →∞ P (N ) d (e -s/N ) exists. Since each converging subsequence yields the same limit, the limit itself exists (see corollary page 59 in [START_REF] Billingsley | Convergence of Probability Measures[END_REF]), i.e. 1 N (M, Z = d) d → M • 1 {Z =d} , as N → ∞, with the limiting vector as in Eq. ( 24).

This concludes the proof of Theorem 6.1.

Conclusion and future work

We first studied the workload for a queue with modulated arrivals, service requirements and service capacity, and derived that the scaled workload converges to an exponentially distributed random variable in heavy traffic. The workload results obtained are valid for any service distribution and for any service discipline which does not depend on the environment. We then focussed on the special setting of a multi-class queue under the DPS policy and showed that the joint queue-length distribution for such a system undergoes a state-space collapse in heavy traffic. Under the scaling of (1-ρ ∞ ), the vector-valued limiting distribution is independent of the modulating environment and converges in distribution to a one-dimensional random variable times a deterministic vector. In this derivation, the distribution of the scaled workload is a key quantity. With this we extend known results about the DPS queue to a Markov-modulated setting.

Clearly an interesting question for future consideration is whether the state-space collapse for the DPS policy carries over to continuously modulated service requirements, as discussed in Remark 3. Another open question concerns the characterization of the moments of the queue lengths for the modulated DPS queue, outside of heavy traffic. Last but not least, modulating the weights of the DPS would open the possibility of dynamical scheduling based on the environment. The latter would be a study on its own, as already the stability conditions will no longer be independent of the weights of the DPS policy.

Acknowledgments The authors would like to thank Urtzi Ayesta, Joke Blom and Michel Mandjes for helpful discussions. This research was partially supported by the SMI Program of INP Toulouse.

Appendix: Proof of Theorem 3.3

The proof Theorem 3.3 is based on Theorem 4 in [START_REF] Dimitrov | Single-server queueing system with Markov-modulated arrivals and service times[END_REF], which can be adapted to our model as follows:

We start with notation and some preliminaries. c D ) and p 0 = ( p 0,1 , . . . , p 0,D ). Furthermore H1 and H2 are the diagonal matrices corresponding to the moments h d1 and h d2 , respectively, for d = 1, . . . , D. Recall Eq. ( 3),

Let = diag(λ 1 , . . . , λ D ), H (s) = diag(1 -h 1 (s), . . . , 1 -h D (s)), C = diag(c 1 , . . . ,
[Q • a] d = c d -λ d h d1 - c ∞ (1 -ρ ∞ ).
We will now construct a partial inverse of Q to make it easier to find a vector a which solves this equation. Let Q 1 and R be matrices such that

Q 1 = ⎛ ⎜ ⎝ q 22 q 23 • • • q 2D . . . q D2 q D3 • • • q D D ⎞ ⎟ ⎠ , R = ⎛ ⎜ ⎜ ⎜ ⎝ 0 0 • • • 0 0 . . . Q -1 1 0 ⎞ ⎟ ⎟ ⎟ ⎠ .
Then det Q 1 = 0 and, due to Q being a generator (for more details, see [START_REF] Dimitrov | Single-server queueing system with Markov-modulated arrivals and service times[END_REF]), we have

Q R = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 -π 2 π 1 -π 3 π 1 • • • -π D π 1 0 1 0 • • • 0 0 0 1 0 • • • . . . . . . . . . . . . . . . 0 0 • • • • • • 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .
It follows that for any vector x, it holds that

x Q R = x - x 1 π 1 π . ( 33
)
Then it can be verified with straightforward calculations that a = (a 1 , . . . ,

a D ) = R(c 1 -λ 1 h 11 -c ∞ (1-ρ ∞ ), . . . , c D -λ D h D1 -c ∞ (1-ρ ∞ )) T = R[C -H1 ]e -c ∞ (1 -ρ ∞ )Re = R[C -H1 ]e -r ( 34
)
is a possible solution vector, with r := c ∞ (1 -ρ ∞ )Re.

Define the vector ϕ = (ϕ 1 (s), . . . , ϕ D (s)) and write Eq. ( 5) in matrix-vector terms:

ϕ(s)Q = ϕ(s)[ H (s)

-sC] + s p 0 C. [START_REF] Zorine | On ergodicity conditions in a polling model with Markov modulated input and statedependent routing[END_REF] Observe that, according to Eq. ( 7),

p 0 C e = c ∞ (1 -ρ ∞ ).
Now multiply from the right both sides of the new vector-matrix equation, Eq. ( 35), with a D-dimensional vector of 1's, e, to obtain

ϕ(s)[ H (s) -sC]e + sc ∞ (1 -ρ ∞ ) = 0. ( 36
)
Multiply from the right both sides of Eq. ( 35) with the matrix R to get Rearranging Eq. (40) yields 1)

ϕ 1 (s/N ) = π 1 B 1 (s/N ) -(B 2 (s/N ) + B 3 (s/N )) = π 1 c ∞ s/N 2 + o(N -2 ) (s/N ) 2 D d=1 π d [λ d h d2 /2 + a d (λ d h d1 -c d )] + c ∞ s/N 2 + o(N -2 ) = π 1 1 + o(
1 + c -1 ∞ d π d [λ d h d2 /2 + a d (λ d h d1 -c d )] s + o(1)
.

Let M be the desired mean stated in Theorem 3.3, that is,

M := c -1 ∞ d π d λd h d2 /2 + a d ( λd h d1 -c d ) .
Then, taking the heavy-traffic limit,

Remark 1 (

 1 Modulated service requirements rewritten to classes) Any multi-class system where the service requirement distribution is determined by the modulating process at a customer's arrival can be written as a multi-class model with nonmodulated, only class-dependent, service rates μ k , as illustrated below: While in state d of the environment, class-k customers arrive with rate λ k,d and have exponential service requirement with mean 1/μ k,d and weight g k . Such customers we refer to as being of class (k, d) and count with M k,d ; hence, we need to keep track of K • D different customer "classes". Arrivals to class (k, d) are inactive while not in state d.

d z m , where z m 1 1 .

 1 . . z m K K =: z m and (m 1 , . . . , m K ) ≥ (0, . . . , 0), i.e. m ≥ 0. Then the overall generating function for the queue length is P(z

2 + 2 = -(s/N ) 2 D d=1 π d λ d h d2 / 2 - 2 H=

 22222 ϕ(s)Q R = ϕ(s)[ H (s) -sC]R + s p 0 C R.Rewrite this equation by using the property of Eq. (33) to obtainϕ(s) = ϕ 1 (s) π 1 π + ϕ(s)[ H (s) -sC]R + s p 0 C R. (37)Iterate Eq. (37) with itself by inserting ϕ(s) into the right-hand side of the equation to obtain, after some algebraic transformations,ϕ(s) = ϕ 1 (s) π 1 π [I + G(s)R] + y(s), (38)with G(s) := H (s) -sC andy(s) := ϕ(s)[G(s)R] 2 + s p 0 C R[G(s)R + I ].(39)Substitute Eq. (38) into Eq. (36) to obtain an expression for ϕ 1 (s):0 = ϕ 1 (s) π 1 π[I + G(s)R] + y(s) • G(s)e + sc ∞ (1 -ρ ∞ ) = ϕ 1 (s) π 1 [π G(s)e + π G(s)RG(s)e] + y(s)G(s)e + sc ∞ (1 -ρ ∞ ) = ϕ 1 (s) π 1 [B 2 (s) + B 3 (s)] + B 1 (s), (40) with B 1 (s) = y(s)G(s)e + sc ∞ (1 -ρ ∞ ), B 2 (s) = πG(s)e and B 3 (s) = πG(s)RG(s)e.The next step is to insert the scaling s → s/N for each term. Recall that using the heavy-traffic parametrization introduced in Sect. 2,wehave (1 -ρ ∞ ) = 1/N.Now observe that, as N →∞, H (s/N ) , and since |ϕ(s/N )| ≤ 1 and p p(N ) 0 C R[G(s/N )R + I ] → 0,as N → ∞. Combining the above we obtainB 1 (s/N ) = y(s/N )G(s/N )e + sc ∞ (1 -ρ ∞ )/N = sc ∞ (1 -ρ ∞ )/N + o(N -2 ) = sc ∞ /N 2 + o(N -2 ).ThenB 2 (s/N ) = π H (s/N ) -sC/N e = π H (s/N ) -H1 s/N e + π H1 s/N -sC/N e = π H (s/N ) -H1 s/N (s/N ) 2 (s/N ) 2 esc ∞ (1 -ρ ∞ )/N = -π H2 2 (s/N ) 2 e + o(N -2 )sc ∞ /N sc ∞ /N 2 + o(N -2 ), since π H1 e = D d=1 π d λ d h d1 = ρ ∞ c ∞ . Furthermore, B 3 (s/N ) = π H (s/N ) -Cs/N R H (s/N ) -Cs/N e = π (s/N ) π (s/N ) 2 H1 -C R H1 -C e + o(N -2 ) = -π (s/N ) 2 H1 -C (a + r) + o(N -2 ) =d π d [(a d + r d )(λ d h d1c d )] (s/N ) 2 + o(N -2 ),due to R C -H1 e = a + r; see Eq.[START_REF] Yechiali | Queueing problems with heterogeneous arrival and service[END_REF]. Under the heavy-traffic scaling,r = c ∞ (1 -ρ ∞ )Re = c ∞ N -1 Re is an o(1) term. Observe that -(B 2 (s/N ) + B 3 (s/N )) = (s/N ) 2 D d=1 π d [λ d h d2 /2 + (a d + o(1))(λ d h d1c d )] + sc ∞ /N 2 + o(N -2 ).

  Ms ,i.e. the LST ϕ(s) converges in distribution to the LST of an exponentially distributed random variable with mean M.

Table 1

 1 A multi-class system where the service requirement distribution is fixed upon arrival can be translated into one with non-modulated service requirements

	Classes	Arrival rates	Serv. rate	Weight

  3, we first derive the service requirement of an arbitrary customer while in state d. If H k (•) is the distribution function of a classk customer's service requirement, then the probability of a class-k customer arriving and requiring service not exceeding x is α k,d H k (x). Summing over k now yields the desired distribution,