Dr Bo Hong

Nauth Panday

Jian Shen

Harry V Wang

Wenping Gong

Anna Soehl

Modeling water exchange between Baltimore Harbor and Chesapeake Bay using artificial tracers: seasonal variations

Keywords: Water exchange, three-dimensional model, Tracer, Age, Baltimore Harbor, Chesapeake Bay

Understanding the dynamics of water exchange between Baltimore Harbor and the Chesapeake Bay is essential when evaluating transport and fate of dissolved substances in both of these systems. Conservative artificial tracers are used in this study to investigate transport processes through a three-dimensional hydrodynamic model (CH3D). The model well reproduced the three-layered circulation pattern in Baltimore Harbor. Several numerical experiments are performed to trace the water mass coming from different sources. The results indicate that both the upper and lower layers of the Harbor are the dominant pathways of transporting dissolved substances from Susquehanna River to the Harbor. Such inward transport is intensified (suppressed) during the high-discharge (low-discharge) period. The upper layer inflow transports water mass with high concentrations of dissolved substances while the inflow from the lower layer transports water mass with low concentrations of dissolved substances. The bottom layer is the dominant pathway for transporting dissolved substances from the lower Bay to the Harbor. Lower river discharge and stronger along-Bay pressure gradient (resulting in stronger landward residual flow in the bottom layer of the Bay) facilitate the bottom intrusion of dissolved substances from lower Bay to the Harbor. Once contaminants are transported into the Harbor, they usually stay for a longer time in the mid-depth of the Harbor than those in other layers due to the threelayer circulation in the Harbor. The time needed for the contaminants being transported out of the Harbor during a typical low-discharge period is about 1 month longer than that needed during a typical highdischarge period. The results, from the environmental perspective, provide new insights for quantitative evaluation on the transport processes of the dissolved biogeochemical substances between Baltimore Harbor and Chesapeake Bay.

Introduction

Baltimore Harbor is a tributary embayment in the upper western portion of Chesapeake Bay draining a highly urbanized and industrialized Patapsco River basin (Fig. 1). Water quality and ecosystem health in the Harbor is susceptible to the industrial, municipal, and stormwater runoff from the surrounding watershed as well as the inflow from the adjacent Bay. Excessive nutrient runoff from the Chesapeake Bay drainage basin is known to contribute to the hypoxia or anoxia events in deepwater (referred to as "dead zone") in the mainstem of the Bay during much of the summer [START_REF] Kemp | Eutrophication of Chesapeake Bay: historical trends and ecological interactions[END_REF]. These deepwater are not only harmful to the aquatic life in the Bay itself, but can also spread into most of the nearby embayments [START_REF] Kuo | Hypoxia and salinity in Virginia Estuaries[END_REF]. Owing to the different source and transport pathways of dissolved substances, the quantitative interpretation of water exchange between Baltimore Harbor and the Bay is necessary for better protecting and utilizing the environmental resources in these areas.

The stratification in upper Chesapeake Bay is maintained by dynamic processes associated with freshwater discharged into the system primarily by the Susquehanna River, resulting in a two-layered circulation [START_REF] Schubel | Responses of upper Chesapeake bay to variations in discharge of the Susquehanna River[END_REF]. Baltimore Harbor has relatively small drainage areas. The combined average daily discharge of the tributaries to the Harbor represents only about 1/315 of the volume of the Harbor [START_REF] Schubel | Responses of upper Chesapeake bay to variations in discharge of the Susquehanna River[END_REF]. The mean tidal range of the Harbor is about 0.35 m and the tidal currents in the Harbor are relatively weak (http://tidesandcurrents.noaa.gov). The estimated maximum tidal excursion is about 2.63 km.

Given the insignificant freshwater discharge and the low tidal energy input to the Harbor, the unique three-layered circulation developed in the Harbor and the stratification are mainly controlled by the density conditions in the adjacent Bay.

Although the three-layered circulation in Baltimore Harbor has been extensively investigated by many previous studies (e.g., [START_REF] Carpenter | The Chesapeake Bay Institute study of the Baltimore Harbor[END_REF][START_REF] Sheng | A three-dimensional mathematical model of coastal estuarine and lake currents using a boundary fitted grid[END_REF][START_REF] Cameron | Estuaries[END_REF][START_REF] Chao | Three layered circulation in reverse estuaries[END_REF], quantitative evaluation of the water exchange and pollutant transport processes between the Harbor and the Bay from the environmental prospective is still very limited. The questions that need to be further addressed include: (i) what amount of dissolved pollutants in the Baltimore Harbor is associated with the inflow from the main Bay? (ii) how long does it take for these substances to be transported from the source area to the area of Transport of water and dissolved substances in estuary is strongly influenced by numerous factors such as river discharge, tides, wind, bottom topography, Coriolis effect, etc. It is usually difficult to trace water mass ( or nutrients, dissolved pollutants, etc. ) coming from different source regions in in-situ observations since hydrodynamic processes usually result in the mixture of them. However, the fate and transport timescales of these biogeochemical substances are essential to both oceanographic research and environmental assessment. Using artificial tracers in numerical modeling is one of the feasible methods that can be used to address the problem.

Previous studies have demonstrated the advantages of using artificial tracers in numerical simulations to determine the fate and transport timescales of dissolved substances since they left the source region. For example, [START_REF] Hirst | Determination of water component age in ocean models: application to the fate of North Atlantic Deep Water[END_REF] examined the penetration and fate of North Atlantic Deep Water in a global ocean model; [START_REF] Delhez | The concept of age in marine modeling II. Concentration distribution function in the English Channel and the North Sea[END_REF] simulated the advectiondispersion of tracers discharged at the nuclear fuel reprocessing plant of Cap de La Hague; [START_REF] Shen | Calculating age and residence time in the tidal York River using threedimensional model experiments[END_REF] used a three-dimensional numerical model to calculate the age distribution for the substance released from the head waters of the York River estuary under different hydrodynamic conditions. In addition, both particle trajectories and passive tracers were used in a three-dimensional circulation model by [START_REF] Döös | Assessment of water exchange between a discharge region and the open sea -A comparison of different methodological concepts[END_REF] to estimate the potential fate and distribution of radio-nuclides released in the coastal region of the Baltic Sea. Through age tracer modeling, [START_REF] Gustafsson | Elucidating the dynamics and mixing agents of a shallow fjord through age tracer modeling[END_REF] quantified the timescales of downward (upward) mixing of surface (bottom) water in a shallow fjord. In these studies, artificial tracers are integrated simultaneously with the hydrodynamic fields and provide powerful tools that can help us to quantify the transport processes of dissolved substances and understand the mechanisms that control their temporal and spatial variations.

For the purpose of quantifying water exchange between Baltimore Harbor and the mainstem of the Chesapeake Bay, we add conservative artificial tracers in an extensively calibrated three- Dimensions, CH3D) was used in this study. CH3D was originally developed by [START_REF] Sheng | A three-dimensional mathematical model of coastal estuarine and lake currents using a boundary fitted grid[END_REF], and subsequently modified extensively by the US Waterways Experiment Station for application to Chesapeake Bay [START_REF] Johnson | User's guide for a three-dimensional numerical hydrodynamic model of Chesapeake Bay[END_REF][START_REF] Johnson | Validation of threedimensional hydrodynamic model of Chesapeake Bay[END_REF]1995;[START_REF] Wang | Application of vertical turbulence closure Schemes in the Chesapeake Bay Circulation Model -A comparative study[END_REF]. The detailed description and theoretical aspects of CH3D can be found in [START_REF] Johnson | User's guide for a three-dimensional numerical hydrodynamic model of Chesapeake Bay[END_REF]. Second order k-ε turbulence closure model (Kundu, 1980) is used to calculate the eddy viscosity and diffusivity coefficient. The background horizontal diffusivity is set to 1×10 -5 m 2 s -1 . A particular feature of CH3D is the solution of transformed equations on a boundary-fitted curvilinear coordinates in the horizontal plane. This allows the model grid to cope with the complicated shoreline configuration and keep a good representation of the deep channel orientation in the Chesapeake Bay and its tributaries, whereas high-resolution grids are required in other models with orthogonal curvilinear coordinates to meet the needs of keeping reasonable orthogonalization and good representation of complex estuarine bathymetry.

CH3D has been successfully used to investigate the hydrodynamic features and transport processes of dissolved substances in the Chesapeake Bay (e.g., [START_REF] Hood | Modeling particles and pelagic organisms in Chesapeake Bay: Convergent features control plankton distributions[END_REF][START_REF] Xu | Assimilating high-resolution salinity data into a model for a partially mixed estuary[END_REF][START_REF] Shen | Determining the age of water and long-term transport timescale of the Chesapeake Bay[END_REF]. It is also extended to include the sediment transport [START_REF] Lin | Toxic modeling in an industrial Harbor -A case study for Baltimore Harbor[END_REF][START_REF] Park | A model study of the estuarine turbidity maximum along the main channel of the Upper Chesapeake Bay[END_REF] and biogeochemical sub-model [START_REF] Xu | Modeling biogeochemical cycles in Chesapeake Bay with a couple physical-biological model[END_REF], in which the physical model provides hydrodynamic background for the sediment and biogeochemical constituents in three dimensions. The model grid of this application (developed by [START_REF] Wang | The development of a water quality model for Baltimore Harbor, Back River, and the adjacent upper Chesapeake Bay. A report to Maryland Department of the Environment[END_REF] is shown in Fig. 2. There are 19 layers in the vertical with a uniform layer thickness of

Methods

Water coming from a specific origin is identified by the passive tracer, which also represents the dissolved substances (contaminants, etc.) in the water. There are two aspects that need to be addressed for quantifying the exchange of water coming from different origins. One is the amount of the dissolved substances, which is represented by the tracer concentration. The other is the transport timescale of the dissolved substances, which is represented by the tracer age (the tracer age is defined as the time elapsed since the tracer under consideration is released into the water body where its age is prescribed to be zero [START_REF] Delhez | Toward a general theory of the age in ocean modeling[END_REF]). The concentration and the age of the tracer are calculated in the numerical model using transport equations [START_REF] Deleersnijder | The concept of age in marine modeling I. Theory and preliminary model results[END_REF]. The tracer concentration equation is

( , ) ( ( , ) ( , )) 0 c t x vc t x K c t x t ∂ + ∇ -∇ = ∂ r r r r (1)
Where ( , ) c t x r is the tracer concentration, v r is the velocity vector, K is the diffusivity tensor, t is the time, and x r is the spatial coordinate. There is no sink of the tracer within the study area.

According to [START_REF] Deleersnijder | The concept of age in marine modeling I. Theory and preliminary model results[END_REF], the evolution of age concentration ( ( , )

t x α r ) is described as ( , ) ( ( , ) ( , )) ( , ) t x v t x K t x c t x t α α α ∂ + ∇ -∇ = ∂ r r r r r (2)
Then the mean age ( ( , ) a t x r ) can be calculated as

( , ) ( , ) ( , ) t x a t x c t x α = r r r (3)
At the releasing site, the tracer concentration is 1 (arbitrary unit) and the tracer age is zero. The tracer age at any other location is representative of the timescale for water or a dissolved

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 8
substance to be transported from its source to this location. The spatial and temporal variations of tracer concentration and mean age (hereafter referred to "water age") will be used in the following to quantify the water exchange under different hydrodynamic conditions.

Numerical Experiments Design

In this study, tracers are released from three different locations in three independent experiments for quantifying water exchange between the Baltimore Harbor and Chesapeake Bay.

Previous studies have indicated that the circulation in the upper Chesapeake Bay is driven primarily by the discharge of Susquehanna River [START_REF] Schubel | Responses of upper Chesapeake bay to variations in discharge of the Susquehanna River[END_REF]. In order to obtain more general results of seasonal variation, we use the climatological daily river runoff data instead of using the data from a single year. In reality, the perturbation/mixing induced by wind forcing is more complicated [START_REF] Maccready | A model study of tide-and windinduced mixing in the Columbia River Estuary and plume[END_REF]. We exclude the wind effect in this study for simplicity and only focus on the effects of river discharge, which is the dominant physical process determining the seasonal variation in the upper Bay. We are aware of the importance of wind mixing as it affects the estuarine stratification. Wind effects are more dominant during synoptic events. The influences of short-term wind events on long-term stratification and long-term transport warrant more observations and studies and will be studied Hydrodynamic equations and passive tracer equations are integrated simultaneously. In order to reach the dynamic equilibrium condition and eliminate transition errors resulting from model spin-up, the model is run repeatedly for 2 years with the annual cycle forcing fields. Only outputs of the last year are used for the following analyses. Note that the spin-up time in Chesapeake Bay is relatively longer than most estuaries, for example, the spin-up time is less than two weeks for Columbia River estuary hindcast modeling [START_REF] Liu | Evaluation of a coastal ocean circulation model for the Columbia River plume in 2004[END_REF]. This is due to the relatively slow adjustment of baroclinic fields in the highly stratified Chesapeake Bay.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 9

Circulation patterns

To understand the long-term variability of water exchange between the Harbor and the Bay, two fundamental processes need to be addressed. The first one is the variation of the Susquehanna River discharge. Susquehanna River has a seasonal flow pattern typical of midlatitude rivers: high discharge in spring produced by snow melting and spring rains followed by low-to-moderate flow throughout most of the remainder of the year [START_REF] Carpenter | The Chesapeake Bay Institute study of the Baltimore Harbor[END_REF][START_REF] Schubel | Responses of upper Chesapeake bay to variations in discharge of the Susquehanna River[END_REF]. This variation greatly impacts the vertical stratification of the main Bay, horizontal salinity gradients between the Bay and Baltimore Harbor, and the amount of dissolved substances transported by river discharge. The second one is the impact of deep water coming from lower Chesapeake Bay, in which the high-salinity water is advected up-Bay mainly through non-tidal gravitational circulation [START_REF] Pritchard | The dynamic structure of a coastal plain estuary[END_REF]. Because of the persistence of stratification and strong two-layered circulation developed in the main Bay, the flow in the bottom layer not only plays a role of conveying saline water from the lower Bay to the upper Bay, but also conveys biogeochemical substances that mix into the deep water during various destratification events.

The Susquehanna River discharge accounts for 87% of the total river runoff in the upper Chesapeake Bay. The 23-year (1985-2007) daily mean river runoff of Susquehanna River (Fig. 4) indicates that the maximum (minimum) river discharge usually appears in March-April (July-October), which can be defined as the typical high-discharge (low-discharge) period. Baltimore

Harbor is famous on its three-layered circulation pattern driven by differences in the vertical salinity structure in the Harbor and in the adjacent Chesapeake Bay [START_REF] Carpenter | The Chesapeake Bay Institute study of the Baltimore Harbor[END_REF][START_REF] Sheng | A three-dimensional mathematical model of coastal estuarine and lake currents using a boundary fitted grid[END_REF][START_REF] Boicourt | A hydrodynamic study of the Baltimore Harbor system[END_REF][START_REF] Chao | Three layered circulation in reverse estuaries[END_REF]. The typical three-layered circulation in Baltimore Harbor is shown in Fig. 5a (after [START_REF] Cameron | Estuaries[END_REF]. The characteristic salinity distribution of the Harbor has the conventional positive seaward gradient in the lower layer, but a negative seaward gradient nearer the surface. Thus, there is flow into the Harbor both at the surface and at the bottom driven by the pressure gradient. These inflows return to the Chesapeake Bay at mid-depth. The mean salinity of Chesapeake Bay is a product of two competing processes, i.e., freshwater influx into the Bay through rivers and exchange with the adjacent ocean [START_REF] Austin | Estimating the mean ocean-bay exchange rate of the Chesapeake Bay[END_REF]. The

Results

The introduction of artificial tracers in a numerical model is a convenient way to quantify the timescale for the transport process of dissolved substances relative to its source region. It also River, the age distribution in Fig. 8 represents the transport timescale needed for the dissolved substances to reach the Harbor. The common feature of the age distribution is that the water age at either the surface or bottom layer of the Harbor is usually lower than that in the mid-depth of the water column. The largest vertical age difference can be about 50 days and appears in winter.

After the high-discharge period starts, such vertical difference decreases to about 20 days. The lowest water age appears in April when Susquehanna River is at its maximum discharge. The age starts to increase from May and reaches its maximum in November.

More detailed structure of the water exchange process can be observed from the time evolution of tracer concentration and water age at selected stations (Fig. 9). Stations C, D, and E are located in the mouth, the exterior, and the interior of the Harbor respectively, and Station B is in the main channel of the Bay (see Fig. 2 for these locations). For the water mass coming from the Susquehanna River, the tracer concentration decreases with depth in the Harbor (Figs. 9a-c). The highest tracer concentration (>80% at the mouth of the Harbor, >60% in the interior of the

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 12
Harbor) appears in April. The time required for the tracer concentration in the lower layer to reach its maximum value (>30%) is delayed approximately a half-month compared with that in the surface layer. When examining the tracer concentration in the main channel of the Bay (Station B, Fig. 9d), we can find such temporal and spatial variations of tracer concentration also occur in the main Bay. The corresponding water age at each station in the Harbor (Figs. 9e-g)

show significant temporal variations. In general, the water age during the high-discharge period is about 2 months longer than that during the low-discharge period. As the inflow is greatly weakened during the low-discharge period (Fig. 6 and Fig. 7f), the bottom layer water age in the inner Harbor is aproximately 1 month larger than that at the surface layer due to reduced gravitation circulation. The largest water age appears at about 6 m (4 m) depth during the highdischarge (low-discharge) period, corresponding to where the mid-depth outflow exists. Such pattern disappears in the main channel of the Bay (Station B, Fig. 9h) where the water age always increases with depth due to the stratification effect (Shen and Lin, 2006). where R is the total tracer concentration and t τ is the turn over time. So the total tracer concentration will be reduced by a factor of e -1 at the time t τ τ

= . The turnover time equals the residence time if the concentration in an embayment decrease exponentially [START_REF] Zimmerman | Mixing and flushing of tidal embayments in the western Dutch Wadden sea, Part I: distribution of salinity and calculation of mixing time scales[END_REF].

The turnover time of the Harbor is 21.7 days (39.5 days) under the high-discharge (lowdischarge) condition.

Discussion

The results presented in the previous section (summarized in Table 1) reveal that tracers coming from different origins have different pathways and different transport timescales under the three-layered circulation condition in the Baltimore Harbor. These reflect the difference of water exchange processes with respect to different source origins.

Dominant Pathway of Water Exchange for Different Water Sources

The dissolved substances that originate from Susquehanna River (lower Bay) mainly occupy the upper (lower) layer of the water column both in the upper Bay and in the Harbor and their concentrations decrease (increase) consistently with depth (Figs. 9 and 11). However, this kind of symmetric distribution does not appear in the vertical age distribution. In the case of releasing tracers from the Susquehanna River (Exp. 1), the water age profile along the navigation channel of the Harbor shows a 'sandwich' structure (high water age in the mid-depth and low water age in both the upper and the lower layers, Fig. 8). For the case of releasing tracers from the lower Bay (Exp. 2), the water age within the Harbor shows a monotonic vertical structure (decreasing consistently with depth, Fig. 10). the transport of water that originates from the lower Bay through the upper layer is still much slower than that in the lower layer, even during a typical high-discharge period with strong surface inflow to the Harbor.

Comparison with theoretical scales

Examining the age distribution for Exp. 2 (Fig. 11), the low age distribution near the bottom does not follow the freshwater discharge, which can be attributed to the influence of the lower Bay through the gravitational circulation. The observed salinity values at Station CB4.4 (see Fig. by about 8 cm s -1 (Fig. 15).

The gravitational circulation in the estuary is influenced by the combined effect of river discharge and bottom saline water intrusion from the ocean. The competition between freshwater inflow from the head and high-salinity water intrusion from the mouth results in the variation of the along-Bay pressure gradient. The [START_REF] Hansen | Gravitational circulation in straits and estuaries[END_REF] solution (HR65 hereafter) for the along-channel velocity provides a classical way to quantify the contribution of along-channel pressure gradient and river discharge to the gravitational circulation. Progress has been made since HR65 and several assumptions in HR65 have been relaxed now (e.g., [START_REF] Geyer | The dynamics of a partially mixed estuary[END_REF][START_REF] Maccready | Toward a Unified Theory of Tidally-Averaged Estuarine Salinity Structure[END_REF][START_REF] Maccready | Estuarine Adjustment[END_REF]. MacCready and Geyer ( 2010) provide an excellent review on the recent advances of the theoretical solution and formulae for calculating adjustment timescales under perturbations. Using the solution technique of MacCready ( 2007), the along-channel velocity ( , , ) u x z t can be divided into depth-averaged (overbar) and deviations from depth- averaged (prime) components so that u u u′ = + . The quantity

r u Q A =
where r Q is the river discharge and A is the cross-sectional area. The quantity u′ has the standard vertical profile as described in MacCready ( 2004):

{ { 1 2 ( ) ( ) ' E a b u uP u P = + Where 2 1 1 3 2 2 z P h   = -    and 2 3 2 1 9 8 z z P h h     = - -         (5)
The magnitude of the exchange flow is ( ) The depth-averaged along-estuary salinity gradient x s is set to be constant in the HR65 solution.

3 48 E x m u g s H K β = , H is the channel depth,
Here we use the solution of MacCready (2004) (Eq. 7) which allows x s to vary along the estuary.

Following MacCready (2007), we assume that

0 m d t K a C U H = , where d C is a drag coefficient M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 18 (= 3 2.6 10 - ×
) and 0 a is a tuning coefficient (0.0325 in MacCready ( 2007), but 0.0227 in this study). Part (a) in Eq. ( 5) represents the circulation contributed by river discharge (RV), and part (b) represents the circulation contributed by pressure gradient (PG). Channel depth and width are prescribed based on the realistic conditions in the upper Chesapeake Bay. In the following, we use ocn s to represent the oceanic salinity, which are based on observations at Station CB4.4 (Fig.

15

.1). Based on our 3D model configuration, we choose river discharge values of 2400 m 3 s -1

and 800 m 3 s -1 , which largely represent the conditions of high-discharge and low-discharge periods, respectively. Using the linear adjustment time theory of MacCready ( 2007), the salt intrusion length can be scaled by

0 2 R ocn s L dx s ∝ ∫
and the freshwater replacement time can be approximated by 1 2 L u . When river discharge is 2400 m 3 s -1 (800 m 3 s -1 ), the freshwater replacement time is about 15.2 (56.0) days, and the salt intrusion length is about 104 km (128 km). The estimated freshwater replacement timescales are quite consistent with the water ages at the mouth calculated by Exp.1 which represents the mean transport timescale of the freshwater from the Susquehanna River. As the freshwater moves downstream, it replaces the old water.

For transport of water from the lower Bay, the transport timescale depends on the exchange flow.

The Knudsen relation provides a suitable estimation of the exchange flow from a given stratification. For that could be supported in the estuary given the density difference between ocean and freshwater.

Following the methodology of MacCready and Geyer (2010) and using the theoretical solution, the graphical presentation of these scaling results is shown in the plane of ( E u c , Φ ) (Fig. 16).

The Knudsen relation is purely diagnostic. It tells us only which exchange flow is consistent with a given stratification. On its own, it cannot predict both stratification and exchange. Meanwhile, our 3D model results are sampled to calculate the dimensionless exchange flow and stratification

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 19
under various real forcing conditions. The 3D model results largely match the theoretical results (Fig. 16). It is indicated that in the upper Chesapeake Bay stronger exchange flow corresponds to stronger stratification. During the high salinity water intrusion events, the horizontal salinity gradients increase and both the vertical stratifications and exchange flow becomes stronger.

Therefore, the age of the water mass coming from thelower Bay is relatively smaller.

Renewal of Baltimore Harbor

Persistent low DO during the summer in the deep-water (dead zone) of the main Bay remains a serious concern for water quality managers. Observations indicate that the hypoxic or anoxic water in the deep channel of the main Bay can enter most of the major tributaries due to gravitational circulation [START_REF] Kuo | Hypoxia and salinity in Virginia Estuaries[END_REF]. This phenomenon has been observed in the Baltimore Harbor, which greatly degrades the ecosystem in the inner Harbor. In the experiment of tracing deep water coming from the lower Bay, we find the upper layer tracer concentration near the head of the Harbor is higher than that in the inner and outer Harbor (figure not shown).

The age distribution (Figs. 7c-d and Fig. 10) suggests that, in the upper layer of the water column, the time needed for the deep water of the lower Bay to reach the head of the Harbor is less than that to the inner and outer Harbor. According to [START_REF] Cameron | Estuaries[END_REF], the characteristic salinity distribution of the Harbor has a positive seaward gradient in the lower layer, but a negative seaward gradient nearer the surface. The resulting stratification is weakest near the head of the Harbor, where vertical mixing is stronger than in other areas of the Harbor.

The increased vertical mixing (or upwelling) near the head of the Harbor may facilitate the upward transport of dissolved substances in that area.

Although the concentrations of riverine pollutants entering the Harbor are higher during the high-discharge period than that during the low-discharge period, the renewal of water in the Harbor is faster during the high-discharge period than that during the low-discharge period.

Thus, the impact of dissolved riverine pollutants on the Harbor is not severe during the highdischarge period as large amounts of substances will be transported out of the Harbor quickly.

The worst situation can occur at the end of the high-discharge period when the concentrations in the Harbor are still high but the water renewal is slowing down. The water exchange between the Baltimore Harbor and Chesapeake Bay is investigated using artificial tracer in conjunction with a calibrated realistic 3D model. The tracer concentration is representative of water or dissolved substances coming from the concerned source region, while the water age can indicate the transport timescale of these dissolved substances. These features are important factors that need to be examined when trying to protect the tidal environmental resources and manage pollutant transport in estuaries.

Conclusions

The results suggest that the temporal and spatial variations of water exchange between the

Harbor and the Bay have obvious seasonal features, which are highly correlated with the Susquehanna River discharge. The hydrodynamic condition during the high-discharge period facilitates the intrusion of dissolved substances coming from Susquehanna River to the Harbor (mainly via the upper layer of the water column). For the case of releasing tracers from Susquehanna River (Exp. 1), the tracer concentrations increase approximately 40%, but the water age decreases approximately 2 months during a high-discharge period compared to the results during a low-discharge period. Both the upper and lower layers are dominant pathways for transporting dissolved substances from the Susquehanna River to the Harbor, but the upper layer tends to transport water masses with higher concentrations of dissolved substances.

During a high-discharge period, the salinity decreases in the upper Bay and the baroclinic pressure between the Bay and the Harbor decreases, which suppresses the saline water intrusion to the Harbor thus decreases the transport time of water coming from the lower Bay (Exp. 2).

Compared to a high-discharge period, the tracer concentration decreases about 20% and the water age increases about 20 days during a low-discharge period. Only the lower layer is the dominant pathway of transporting dissolved substances coming from the lower Bay into the Harbor.

The transport in the main channel is highly related to the dynamic conditions in the lower Bay. The long-term mean transport processes, which mainly depend on the low frequency and mean motion of the water, are essential to evaluate water exchange between the Harbor and the Bay.

This study indicates that the water age and tracer concentration in the Harbor are sensitive to the perturbations of the Susquehanna River discharge. In reality, the perturbation can also be induced by the wind forcing. and affect the water exchange over shorter timescales. The incorporation of the wind effect will be covered in our future work. As the freshwater discharge dominates stratification in the upper Bay, the results presented in this study represents the seasonal density-driven exchange flow between the Harbor and the Bay. The approach used here can be used to study the water exchange processes in other estuaries and harbors. Figure 10. Same as Fig. 8, but for the results of Exp. 2. Figure 11. Same as Fig. 9, but for the results of Exp. 2. Figure 12. Temporal evolution of total tracer concentration (arbitrary unit) in Baltimore Harbor under high-flow (thin line) and low-flow (bold line) condition, respectively. The dashed lines indicate where the total tracer concentration is reduced to a factor of e -1 (i.e., measure of the turnover time of the Harbor waters).

  (iii) how does the transport time-scale change spatially and temporally under different hydrodynamic conditions?

  dimensional hydrodynamic model. The hydrodynamic model and the method of quantifying water exchange are introduced in Section 2. Section 3 simply describes the hydrodynamic background of the concerned area. The diagnostic and quantitative analyses of water exchange between Baltimore Harbor and main Bay are presented in Section 4. Discussion and conclusions are provided in Section 5. dimensional numerical hydrodynamic model (Curvilinear Hydrodynamics in Three

  schemes are adopted for the temperature and salinity fields at the southern open boundary condition in which the observed vertical salinity and temperature profiles at Station CB4.4 (see Fig. 2 for its location) are advected into the model domain for the inflow conditions. Since the observed salinity and temperature data are collected at approximately monthly intervals, the data were first linearly interpolated into daily intervals to obtain the annual cycle. Then the station data were interpolated onto the model grid as open boundary forcings. The radiation boundary condition was used for the velocity field. The velocity at the open boundary was calculated in the model based on the observed surface elevation field and the density profile. The daily averaged surface heat exchange and equilibrium temperature are computed from the corresponding meteorological data (NOAA, 2007). Lin et al. (2004) used this model configuration to simulate hydrodynamic fields covering both very wet and very dry hydrological years. Modeled bottom velocity fields were compared with observations at a mooring station inside the Baltimore Harbor in spring of 2000. The simulated along-Harbor velocity agrees very well with the observations. Park et al. (2008) also used this model configuration to simulate the hydrodynamic field in 1996. By using the two surveys in October followed a peak discharge of 5,267 m 3 s -1 on October 21 and the two surveys in July under the low discharge conditions (data presented in Sanford et al. (2001)), the model performance were validated over a wide range of freshwater discharge values. The spatial distribution characteristics of observed along-Bay salinity profiles were reproduced well by the model. The model-data comparisons at a mooring station indicated that the model could catch the observed vertical structure variations and the vertical mixing scheme worked properly. These model-data comparisons indicated the success of model performance on tidal timescales. Liu et al. (2009) presented a good example of model-data comparisons. They evaluated model skill not only for tidal properties but also for subtidal properties. Since we are focusing on the seasonal variations in this study, further model validation has been conducted to confirm the model performance in modeling seasonal estuarine circulation. Whereas simulations were conducted for multiple years, only the results in 2000 are presented here (similar model performances were obtained in other years). Comparisons between the modeled salinity and CBP monitoring data are shown in Fig. 3a. Model results show consistent seasonal variations with the observations. It is shown that the model can catch the seasonal salinity variations in the upper Chesapeake Bay.

  separately. The three tracer simulation experiments are hereafter referred to as Exp. 1, Exp. 2, and Exp. 3. In Exp. 1, tracers are continuously released from the entrance of Susquehanna River, which are used to track the water masses (contaminants, etc.) coming from Susquehanna River.In Exp. 2, tracers are continuously released from the mid-Bay (along the southern open boundary of the model), which are used to track the water masses coming from the lower Bay. In Exp. 3, tracers are instantaneously released over the entire Baltimore Harbor, which are used to track the dissolved substances over the Harbor. There is no other source of tracer in each experiment.

  Figs. 5b and 5c present the typical three-layered circulation (along the navigation channel of the Harbor, transect I in Fig.2) simulated in this study. The model results show similar pattern with that demonstrated by[START_REF] Cameron | Estuaries[END_REF]: in column, the less saline, less dense surface waters of the Bay tend to overflow the more saline, denser surface waters in the Harbor. Near the bottom, the more saline, denser waters of the Bay tend to flow under the less saline, less dense waters of the Harbor. The mid-depth water flows outward to offset the inward flow.

  modeled surface salinity fields in the upper Chesapeake Bay averaged during the typical highdischarge and low-discharge periods are shown in Fig.6a. Competitions between the river runoff and lower Bay saline water intrusion are clearly reflected. During the high-discharge period, water in the upper Bay is greatly diluted. The intrusion of the freshet into Baltimore Harbor can be clearly discerned from the isohalines. During the low-discharge period, the saline water takes over much of the domain gradually and the surface fresher water intrusion from the main Bay to the Harbor is greatly weakened and mainly confined to the area along the northern bank of the Harbor. The residual velocity profile across the mouth of Baltimore Harbor (facing west, transect II in Fig.2) is presented in Fig.6b. Without the wind forcing, three-layered circulation appears throughout the year, but both the spatial distribution and the strength of the flow vary with time.Uniform surface inward flow usually occurs during the high-discharge period, whereas surface inward flow only occurs near the northern bank of the Harbor during the low-discharge period resulting in the two-layered circulation near the southern bank. The mid-depth outward flow is strengthened when the surface inflow is fully developed during the high-discharge period, and then weakened during the low-discharge period. The inward flow in the bottom layer shows less variation compared with those in the upper layers. This is consistent with the findings of[START_REF] Chao | Three layered circulation in reverse estuaries[END_REF] that indicated the existence of persistent and dominant intrusion into the bottom of the Harbor against all adversities.

  can reveal the spatial distribution of tracer concentrations of dissolved substances coming from a M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 11 specific source. These are usually difficult to measure in observations. Water exchange between Baltimore Harbor and Chesapeake Bay are investigated below from three different aspects. 4.1 Water Mass Coming from Susquehanna River (Exp. 1) Exp. 1 tracks the dissolved substances coming from Susquehanna River. The selected snapshots (daily mean results sampled at Day 91 and Day 255, respectively) of surface water age distribution are shown in Figs. 7a and 7b, which correspond to the very high-discharge and very low-discharge conditions (see Fig. 4), respectively. The water age at Day 91 is much lower than that at Day 255, indicating that the transport time of the water mass coming from Susquehanna River during a high-discharge period is much less than that during a low-discharge period. The water age difference between Days 91 and 255 is approximately 70 (40) days at the mouth (head) of the Harbor. The corresponding surface velocity field (Figs. 7e and 7f) shows that the flow at Day 91 is much stronger than that at Day 255. There is strong and consistent surface inflow in the Harbor at Day 91, whereas weak and divergent flow appears in the Harbor at Day 255. These demonstrate that Susquehanna River discharge has substantial influence on water exchange between the Harbor and the main Bay.

Fig. 8

 8 Fig. 8 presents the monthly mean water age profiles along the navigation channel in Baltimore Harbor. Since the tracers are continuously released at the entrance of Susquehanna

4. 2

 2 Fig.11presents the depth-time profiles of tracer concentration and water age at selected stations. Tracer concentrations in the Harbor (Figs.11a-c) increase with depth throughout the year. The lowest tracer concentration appears during the high-discharge period (at both the upper and the lower layers of the water column), and then increases gradually as the low-discharge period begins. The tracer concentration during the low-discharge period is usually higher than that during the high-discharge period by 20-30%, and the vertical difference of tracer concentration during the high-discharge period is higher than that during the low-discharge period. Similar patterns can be seen in the main Bay (Fig.11d). The time evolutions of vertical age profiles at selected stations are shown in Figs.11e-h. The highest vertical age variations appear at around 6 m, where the haloclines fluctuate with the dynamics. The bottom layer water age is almost the same between the Harbor and the Bay. But in the upper layer, the water age in the Harbor is about 30 days larger than that in the main Bay.

  Fig. 13 shows the evolution of tracer concentration in the Harbor under the typical high-discharge condition (a similar pattern develops under the low-discharge condition). It can be shown that dissolved substances in the Harbor usually remain for a longer time in the mid-depth of the Harbor (and upper layer of the inner Harbor) than those in other areas. Thus when organic contaminants are discharged into the Harbor from the drainage basin, the middepth of the Harbor (and the upper layer of the inner Harbor) usually sustains water quality stress over a longer period than other areas. Examining the conditions at Stations f, g, and h can provide additional insight about the influence of Baltimore Harbor outflow on water quality in the main Bay. Stations f and h lie in the south and north of the Harbor, respectively, while Station g lies near the mouth of the Harbor (see Fig. 2). These stations are used to represent different areas in the Bay. The depth-averaged tracer concentrations at each station are shown in Figs. 14a and 14b. Under the high-flow condition, the tracer concentrations at Station g are zero at the beginning and start to increase after the first week. The peak tracer concentration appears at around Day 27, at which point it starts to decrease gradually with small perturbations. Similar temporal variations can be observed at other stations, with peaks occurring at different times at each station (Day 35 and Day 45 for Stations f and h, respectively). The concentration distribution shows the age density function. The calculated mean age values suggest that contaminants in Baltimore Harbor can be transported to the selected stations f, g, and h at around 66, 58, and 73 days, respectively, after being released in the Harbor. The mean age under the low-flow condition is about 1 month longer than that under the high-flow condition (around 97, 82, and 105 days for Station f, g, and h, respectively). The appearance of peak tracer concentration under the low-flow condition is that under the high-flow condition. The time evolution of the vertical tracer concentration profile helps to reveal the distribution of tracer concentration for the entire water column at each station. As an example, the results under the low-flow condition (a similar pattern appears in the high-flow condition) are shown in Fig. 14c. These results suggest that, after the dissolved substances are transported out of the Harbor, they show different vertical occupation patterns in different areas of the main Bay. At the station south (north) of the Harbor, the tracer concentration decreases (increases) with depth. For the station close to the mouth of the Harbor, the highest tracer concentration occurs at around 6 m depth. The difference in the vertical distribution of tracer concentration at different locations suggests that the contaminants in the Harbor may influence different habitats in the mainstem of the Bay.

  The results in Exp. 1 indicate that both the upper and the lower layer inflow of the Harbor are the dominant pathways of transporting dissolved substances from Susquehanna River to the Harbor. The upper layer inflow, associated with the high concentration of dissolved substances, dominates the mass transport during a typical high-discharge period. But through the bottom inflow becomes dominant during a typical lowdischarge period. Although the concentrations of dissolved substances associated with the bottom inflow is less than those in the surface inflow, the transportation by bottom inflow is faster than that by the surface inflow during a typical low-discharge period. In Exp. 2, on the other hand, only the inflow from the lower layer is the dominant pathway for transporting dissolved substances from the lower Bay into the Harbor. Because of the persistence of stratification in the main Bay, which prevents the vertical mixing of bottom water to the surface,

2

  for its location) is presented in Fig. 15a, which were used as the model open boundary conditions in this study. The temporal vertical salinity distribution shows that there are several high-salinity water intrusion events (i.e., salinity > 20 ppt) in the bottom layer of the Bay. The depth-time profiles of water age in Exp. 2 are shown in Figs. 15b and 15c, which are sampled at Stations B and C, respectively. It can be seen that the low water age events (and the high tracer concentration events, (Fig. 11)) are usually associated with the high-salinity water intrusion events. Although Station CB4.4 is miles away from Stations B and C, the bottom layer salinity at Stations B and C (figures not shown) have consistent temporal variations with salinity at Station CB4.4. In order to further examine the underlying transport mechanisms, the depth-time profile of modeled along-Bay velocity at Station V (see Fig. 2) is presented in Fig. 15d. A 15-day lowpass filter was applied to the velocity data before the plotting. Thees results clearly show the two-layered gravitational circulation development in the Bay and the bottom residual flow variations are highly in phase with the bottom salinity field. Velocity profiles at other stations in the main channel of the Bay show similar temporal variations to those at Station V. Those results indicate that the higher bottom salinity field usually corresponds to the stronger landward flow resulting in a decrease of water age in the bottom layer (Figs. 15b-c). The gravitational in changing the transport processes. When the return flow is stronger, tracer transportation is faster, which results in smaller water age. In general, the water age near the bottom of the main Bay decreases about 20 days when the landward residual velocity increases

  g is the acceleration of gravity, and m K is the vertical eddy viscosity. For simplicity, the density anomaly in the solution is assumed to depend only on salinity (

  vertical eddy diffusivity, which has the same definition as Eq. 3.9 in MacCready (2007)). A dimensionless expression for the exchange flow can speed of the fastest internal wave

A

  high-salinity water intrusion event in the bottom layer of the Bay usually corresponds to a low water age event near the bottom both in the Bay and in the Harbor. The transport time depends on the strength of the exchange flow. The water age decreases by approximately 20 days when the landward residual velocity in the main Bay increases by approximately 8 cm/s. Diagnostic analyses suggest that the low river discharge and strong along-Bay pressure gradient facilitate the up-Bay transport of dissolved substance coming from the lower Bay. for the contaminants in the Harbor to be transported to the main Bay during a typical low-discharge period is about 1 month longer than during a typical highdischarge period. Because of the development of the three-layered circulation, dissolved substances usually remain the longest in the mid-depth of the Harbor. The maximum tracer concentration occurs at around 6 m depth in the area close to the mouth of the Harbor. As soon as the dissolved substances are transported out of the Harbor, water mass follows two different pathways. When being transported toward the downstream of the Bay, high tracer concentrations occur near the surface. Conversely, when being transported toward the upstream of the Bay, high tracer concentrations occur near the bottom.

Figure 1 .

 1 Figure 1. Bathymetry of upper Chesapeake Bay. The color gradient represents water depth in meters.Figure 2. Model domain and boundary-fitted computational grids (101 X 153). The red bold lines denote transects along the navigational channel of the Harbor (transect I) and across the mouth of the Harbor (transect II) respectively. The star symbols represent stations where time series analyses are conducted. The rectangle (Station V) in the main channel of the Bay is the sample site of modeled velocity profile. The triangles represent CBP monitoring stations in the Bay.

Figure 3 .

 3 Figure 3. Model-data comparisons. (a) Salinity results. The subscripts s and b represent surface salinity and bottom salinity, respectively ; (b) low-pass filtered water elevation. The observations (black line) are missing in some time intervals. The correlation coefficients (CC) in (b) are at the 95% confidence level.Figure 4. Long-term mean (1985-2007) daily freshwater discharge (x10 3 m 3 s -1 ) from Susquehanna River. Figure 5. (a) Typical longitudinal-vertical section of the salinity (ppt) distribution along the navigational channel of the Baltimore Harbor, with the Chesapeake Bay lies in the right entrance of the figure. Arrows indicate the schematic depiction of the three-layered circulation pattern in the Harbor (after Cameron and Pritchard, 1963). (b) Modeled salinity (ppt) profile along transect I. (c) Modeled velocity (cm s -1 ) profile along section I. The positive (negative) values represent the outward (inward) flow. The abscissa is consistent with model grid coordinates in Figure 2.

Figure 4 .Figure 6 .

 46 Figure6. Surface salinity (ppt) in the upper Chesapeake Bay averaged during (a) the typical highdischarge period and (b) the typical low-discharge period. Zonal velocity (cm s -1 ) profiles across the mouth of Baltimore Harbor (transect II in Fig.2) averaged during (c) the typical high-discharge period and (d) the typical low-discharge period. The positive (negative) values represent the outward (inward) flow. The contour interval is 1 ppt for salinity field and 0.5 cm s -1 for the velocity field. Figure7. Snapshot of surface water age (day) (Exp. 1 (a-b) and Exp. 2 (c-d)) and velocity (cm s -1 ) (e-f) at Day 91 (left) and Day 255 (right). The contour interval is 4 days. Figure8. Monthly mean vertical profile of water age (day) resulting from Exp.1 along the navigation channel of Baltimore Harbor (transect I in Fig.2). The contour interval is 5 days. The abscissa is consistent with the model abscissa in Figure2. Figure9. Temporal vertical tracer concentration (arbitrary unit) (left, a-d) and water age (day) (right, e-h) profiles resulting from Exp. 1 at Stations C, D, E and B (from top to bottom), respectively. The contour interval for the tracer concentration and water age are 0.1 and 10 days, respectively. Figure10. Same as Fig.8, but for the results of Exp. 2. Figure11. Same as Fig.9, but for the results of Exp. 2. Figure12. Temporal evolution of total tracer concentration (arbitrary unit) in Baltimore Harbor under high-flow (thin line) and low-flow (bold line) condition, respectively. The dashed lines indicate where the total tracer concentration is reduced to a factor of e -1 (i.e., measure of the turnover time of the Harbor waters). Figure13. Evolution of tracer concentration (arbitrary unit) in the case of releasing tracers instantly in Baltimore Harbor under typical high-discharge condition (Exp. 3). Chesapeake Bay lies to the right of the figure. The abscissa is consistent with the model abscissa in Figure2. Figure 14. Time series of depth-averaged tracer concentration under (a) typical high-discharge condition and (b) typical low-discharge condition. The temporal vertical profiles of tracer concentration under typical low-discharge condition are shown in (c). In each row, the results at Stations f, g, and h are presented from the left to the right. Figure 15. The temporal vertical profile of (a) observed salinity (ppt) at Station CB4.4; (b) water age (day) at Station C (Exp. 2); (c) water age at Station B (Exp. 2); and (d) residual velocity (cm s -1 ) at Station V. The negative (positive) velocity represents the seaward (landward) flow.Figure16. The plane (uE/c, Φ) presentation of all possible estuary states. Φ and uE/c denotes the dimensionless stratification and exchange flow, respectively. The solid line represents the results calculated using theoretical solution, the triangle represents the results calculated using the 3D model.

Figure 7 .

 7 Figure6. Surface salinity (ppt) in the upper Chesapeake Bay averaged during (a) the typical highdischarge period and (b) the typical low-discharge period. Zonal velocity (cm s -1 ) profiles across the mouth of Baltimore Harbor (transect II in Fig.2) averaged during (c) the typical high-discharge period and (d) the typical low-discharge period. The positive (negative) values represent the outward (inward) flow. The contour interval is 1 ppt for salinity field and 0.5 cm s -1 for the velocity field. Figure7. Snapshot of surface water age (day) (Exp. 1 (a-b) and Exp. 2 (c-d)) and velocity (cm s -1 ) (e-f) at Day 91 (left) and Day 255 (right). The contour interval is 4 days. Figure8. Monthly mean vertical profile of water age (day) resulting from Exp.1 along the navigation channel of Baltimore Harbor (transect I in Fig.2). The contour interval is 5 days. The abscissa is consistent with the model abscissa in Figure2. Figure9. Temporal vertical tracer concentration (arbitrary unit) (left, a-d) and water age (day) (right, e-h) profiles resulting from Exp. 1 at Stations C, D, E and B (from top to bottom), respectively. The contour interval for the tracer concentration and water age are 0.1 and 10 days, respectively. Figure10. Same as Fig.8, but for the results of Exp. 2. Figure11. Same as Fig.9, but for the results of Exp. 2. Figure12. Temporal evolution of total tracer concentration (arbitrary unit) in Baltimore Harbor under high-flow (thin line) and low-flow (bold line) condition, respectively. The dashed lines indicate where the total tracer concentration is reduced to a factor of e -1 (i.e., measure of the turnover time of the Harbor waters). Figure13. Evolution of tracer concentration (arbitrary unit) in the case of releasing tracers instantly in Baltimore Harbor under typical high-discharge condition (Exp. 3). Chesapeake Bay lies to the right of the figure. The abscissa is consistent with the model abscissa in Figure2. Figure 14. Time series of depth-averaged tracer concentration under (a) typical high-discharge condition and (b) typical low-discharge condition. The temporal vertical profiles of tracer concentration under typical low-discharge condition are shown in (c). In each row, the results at Stations f, g, and h are presented from the left to the right. Figure 15. The temporal vertical profile of (a) observed salinity (ppt) at Station CB4.4; (b) water age (day) at Station C (Exp. 2); (c) water age at Station B (Exp. 2); and (d) residual velocity (cm s -1 ) at Station V. The negative (positive) velocity represents the seaward (landward) flow.Figure16. The plane (uE/c, Φ) presentation of all possible estuary states. Φ and uE/c denotes the dimensionless stratification and exchange flow, respectively. The solid line represents the results calculated using theoretical solution, the triangle represents the results calculated using the 3D model.

Figure 13 .

 13 Figure6. Surface salinity (ppt) in the upper Chesapeake Bay averaged during (a) the typical highdischarge period and (b) the typical low-discharge period. Zonal velocity (cm s -1 ) profiles across the mouth of Baltimore Harbor (transect II in Fig.2) averaged during (c) the typical high-discharge period and (d) the typical low-discharge period. The positive (negative) values represent the outward (inward) flow. The contour interval is 1 ppt for salinity field and 0.5 cm s -1 for the velocity field. Figure7. Snapshot of surface water age (day) (Exp. 1 (a-b) and Exp. 2 (c-d)) and velocity (cm s -1 ) (e-f) at Day 91 (left) and Day 255 (right). The contour interval is 4 days. Figure8. Monthly mean vertical profile of water age (day) resulting from Exp.1 along the navigation channel of Baltimore Harbor (transect I in Fig.2). The contour interval is 5 days. The abscissa is consistent with the model abscissa in Figure2. Figure9. Temporal vertical tracer concentration (arbitrary unit) (left, a-d) and water age (day) (right, e-h) profiles resulting from Exp. 1 at Stations C, D, E and B (from top to bottom), respectively. The contour interval for the tracer concentration and water age are 0.1 and 10 days, respectively. Figure10. Same as Fig.8, but for the results of Exp. 2. Figure11. Same as Fig.9, but for the results of Exp. 2. Figure12. Temporal evolution of total tracer concentration (arbitrary unit) in Baltimore Harbor under high-flow (thin line) and low-flow (bold line) condition, respectively. The dashed lines indicate where the total tracer concentration is reduced to a factor of e -1 (i.e., measure of the turnover time of the Harbor waters). Figure13. Evolution of tracer concentration (arbitrary unit) in the case of releasing tracers instantly in Baltimore Harbor under typical high-discharge condition (Exp. 3). Chesapeake Bay lies to the right of the figure. The abscissa is consistent with the model abscissa in Figure2. Figure 14. Time series of depth-averaged tracer concentration under (a) typical high-discharge condition and (b) typical low-discharge condition. The temporal vertical profiles of tracer concentration under typical low-discharge condition are shown in (c). In each row, the results at Stations f, g, and h are presented from the left to the right. Figure 15. The temporal vertical profile of (a) observed salinity (ppt) at Station CB4.4; (b) water age (day) at Station C (Exp. 2); (c) water age at Station B (Exp. 2); and (d) residual velocity (cm s -1 ) at Station V. The negative (positive) velocity represents the seaward (landward) flow.Figure16. The plane (uE/c, Φ) presentation of all possible estuary states. Φ and uE/c denotes the dimensionless stratification and exchange flow, respectively. The solid line represents the results calculated using theoretical solution, the triangle represents the results calculated using the 3D model.

Figure 1 .

 1 Figure 1. Bathymetry of upper Chesapeake Bay. The color gradient represents water depth in meters.

Figure 2 .

 2 Figure 2. Model domain and boundary-fitted computational grids (101 X 153). The red bold lines denote transect along the navigational channel of the Harbor (transect I) and across the mouth of the Harbor (transect II) respectively. The star symbols represent stations where time series analyses are conducted. The rectangle (Station V) in the main channel of the Bay is the sample site of modeled velocity profile. The triangles represent CBP monitoring stations in the Bay.

Figure 4 .Figure 5 Figure 6 .Figure 7 .

 4567 Figure 3. Model-data comparisons. (a) Salinity results. The subscripts s and b represent surface salinity and bottom salinity, respectively ; (b) low-pass filtered water elevation. The observations (black line) are missing in some time intervals. The correlation coefficients (CC) in (b) are at the 95% confidence level.

  1.52 m, except for the top layer whose thickness fluctuates with water level. The model is

	initialized	by	the	Chesapeake	Bay	Program	(CBP)	Monitoring	data
	(http://www.chesapeakebay.net/dataandtools.aspx/).	Daily	fresh	water	inflow
	(http://waterdata.usgs.gov/nwis/dv/) with zero salinity and time-varying temperature was
	prescribed for the major tributaries, which include Susquehanna River, C&D Canal, Chester
	River, Choptank River, and Patapsco River (see Fig. 1 for the locations of these tributaries). At

the southern open boundary the tidal elevations are prescribed by hourly data observed at

Table Table 1

 Table . A summary of the tracer age results in the Baltimore Harbor.

	Description	Exp. 1	Exp. 2		Exp. 3
	Vertical age distribution	Low age occurs at both	Decreases	consistently	Larger in the mid-depth
	pattern	surface and bottom, and	with depth		than those in other areas

Acknowledgements

The funding of the project is provided by the Maryland Department of Environment. This is the contribution number #xxx of Virginia Institute of Marine Sciences, College of William and Mary.

This work was performed using computational facilities at the College of William and Mary which were enabled by grants from Sun Microsystems, the National Science Foundation, and Virginia's Commonwealth Technology Research Fund.


Henri Gouin
email: henri.gouin@univ-cezanne.fr

Liquid-solid interaction at nanoscale and its application in vegetal biology

Keywords: nanofilms, disjoining pressure, cohesion-tension theory, interface motions, Navier length, ascent of sap

Liquid-solid interaction at nanoscale and its application in vegetal biology

Introduction

The model we develop makes it possible to investigate the behavior of the fluids in the nanofilms, and its applications extend to life sciences. A particularly interesting example concerns vegetable biology: the rise and the motion of sap in the highest trees focus many polemics and debates between biologists. Many of them regard the approach known as the cohesion-tension theory (CTT) proposed at the end of the nineteenth century by Dixon and Joly as the only valid one [START_REF] Dixon | On the ascent of sap[END_REF].

As an obvious fact, Flindt reports huge trees as eucalyptus and giant sequoias of more than 130 meters [START_REF] Flindt | Amazing Numbers in Biology[END_REF], but the biophysical determination of maximum size to which trees can grow is not well understood and calculated. The main problem with the understanding of water transport is why the sap is able to irrigate up very high levels in tall trees. The crude sap contains diluted salts but its physical properties are roughly comparable with the water ones. Hydrodynamics, capillarity and osmotic pressure create a ascent of sap of only few tens of meters [START_REF] Zimmermann | Xylem structure and the ascent of sap[END_REF]. To explain the sap ascent phenomenon, Dixon and Joly proposed a cohesion-tension model, followed by a quantitative attempt [START_REF] Van Der Honert | Water transport in plants as a catenary process[END_REF]: liquids are assumed to be subjected to tensions generating negative pressures compensating gravity effects. As pointed out in [START_REF] Zimmermann | Water ascent in tall trees: does evolution of land plants rely on a highly metastable state? Tansley review[END_REF], a turning-point in the confidence of the opponents to the cohension-tension theory for the sap ascent was the experiment which demonstrated that tall trees survive by overlapping double saw-cuts made through the cross-sectional area of the trunk to sever all xylem elements [START_REF] Preston | Movement of water in higher plants[END_REF]. This result confirmed by several authors does not seem in agreement with the possibility of strong negative pressures in microtubes [START_REF] Mackay | The effects of transverse cuts through the stems of transpiring woody plants on water transport and stress in the leaves[END_REF][START_REF] Benkert | Direct measurements of the pressure and flow in the xylem vessels of nicotiana tabacum and their dependence on flow resistance and transpiration rate[END_REF]. Using a xylem pressure probe, the apparatus does not measure any water tension in many circumstances: xylem tension exceeding 0.6 Mpa seems not to be observed and in normal state most vessels may be embolized at a level corresponding about sixty meter height. Moreover, gas-vapor transportation in xylem tubes seems to appear at the top of high trees ( [START_REF] Zimmermann | Water ascent in tall trees: does evolution of land plants rely on a highly metastable state? Tansley review[END_REF] and its references herein). As comments and questions, M.H. Zimmerman wrote in 1983 [START_REF] Zimmermann | Xylem structure and the ascent of sap[END_REF]: "We don't yet fully understand all aspects of xylem-water supply to leaves and have here a wide-open filed of potential very interesting future research. The heartwood is referred to as a wet wood. It may contain liquid under positive pressure while in the sapwood the transpiration stream moves along a gradient of negative pressures. Why is the water of the central wet core not drawn into the sapwood? Free water, i.e. water in tracheids, decreases in successively older layers of wood as the number of embolized tracheids increases. The heartwood is relatively dry i.e. most tracheids are embolized. It is rather ironic that a wound in the wet wood area, which bleeds liquid for a long period of time, thus appears to have the transpiration stream as a source of water, in spite of the fact that the pressure of the transpiration stream is negative most of the time! It should be quite clear by now that a drop in xylem pressure below a critical level causes cavitations and normally puts the xylem out of function permanently. The cause of such a pressure drop can be either a failing of water to the xylem by the roots, or excessive demand by transpiration." Many proponents of the CTT wrote a letter [START_REF] Angeles | The cohension-tension theory[END_REF] to protest against the recent review [START_REF] Zimmermann | Water ascent in tall trees: does evolution of land plants rely on a highly metastable state? Tansley review[END_REF]. They said that "the CTT is widely supported by biological scientists as the only theory consistent with the preponderance of data on water transport of plant". Nonetheless, the problem of possible cavitation in trees remains. Such liquids are strongly metastable and can generate cavitations causing embolisms in xylem tubes made of dead cells [START_REF] Tyree | The cohesion-tension theory of sap ascent: current contreversies[END_REF]. For example, it is interesting to note that in xylem tube -where diameters range between 50 and 400 µm -the crude sap has a surface tension γ lv lower than the surface tension of pure water which is 72 cgs at 20 • Celsius. If we consider a microscopic gas-vapor bubble inside the crude sap with diameter 2 R, the difference between the gas-vapor pressure P vapor and the liquid sap pressure P liquid can be expressed by the Laplace formula: P vapor -P liquid = 2 γ lv /R. But P vapor being positive, unstable bubbles must appear when R ≥ -2 γ lv /P liquid . For a negative pressure P liquid ≤ -0.6 MPa corresponding to more than sixty meters height, we get R ≥ 0.24 µm. In such a case, dynamical bubbles spontaneously appear from germs naturally existing in a crude liquid and cavitation makes the tubes embolized. Consequently, without any biological known process it is difficult to be convinced that xylem tubes are not embolized when they are filled enough with sap up to altitude significantly more important than one hundred meters corresponding to the highest trees.

Our understanding of the ascent and the motion of sap in very high trees differs from the CTT : at a higher level than a few tens of meters -corresponding to the pulling of water by capillary and osmotic pressure -we assume that xylem microtubes are embolized. In addition, we also assume that a thin liquid film -with a thickness of a few nanometers [START_REF] Derjaguin | Surfaces Forces[END_REF][START_REF] Israelachvili | Intermolecular Forces[END_REF] -wets xylem walls up to the top of the tree. At this scale, long range molecular forces stratify liquids and the ratio between tube diameter and sap film thickness allows us to consider tube walls as plane surfaces. In Section 2, using the calculations presented in [START_REF] Gouin | A new approach for the limit to tree height using a liquid nanolayer model[END_REF][START_REF] Gouin | Liquid nanofilms. A mechanical model for the disjoining pressure[END_REF], we reconsider the analytic expression in density-functional theory for a thin heterogeneous liquid film which takes account of the power-law tail behavior dominant in a thin liquid film in contact with a solid [START_REF] Gouin | Energy of interaction between solid surfaces and liquids[END_REF]. The effects of the vapor bulk bordering the liquid film are simply expressed with an other density-functional located on a mathematical surface. With such a functional, we obtained the equations of equilibrium, motion and boundary conditions [START_REF] Gouin | Boundary conditions for a capillary fluid in contact with a wall[END_REF] for a thin vertical liquid film wetting a vertical solid wall and we computed the liquid layer thickness as a function of the film level; these previous results can be extended to mixtures of fluid and perfect gas [START_REF] Gouin | A new approach for the limit to tree height using a liquid nanolayer model[END_REF]. Then, the so-called disjoining pressure of thin liquid layers yields a natural tool for very thin films [START_REF] Derjaguin | Surfaces Forces[END_REF]. The minimal thickness for which a stable wetting film wets a solid wall is associated with the pancake layer when the film is bordering the dry solid wall and corresponds to the maximal altitude [START_REF] Derjaguin | Surfaces Forces[END_REF][START_REF] De Gennes | Wetting: statics and dynamics[END_REF][START_REF] De Gennes | Capillarity and Wetting Phenomena: Drops, Bubbles, Pearls, Waves[END_REF]. The normal stress vector acting on the wall remains constant through the liquid layer and corresponds to the gas-vapor bulk pressure which is currently the atmospheric pressure and consequently, no negative pressure appears in the liquid layer. At the top of very high trees, the thickness of the sap layer is of a few number of nanometers. The negative pressure is only present for the liquid bulk in micropores. Numerical calculations associated with physical values for water yield the maximal film altitude for a wood material corresponding to a good order of the height of the tallest trees.

In Section 3, we consider the flow of sap at high levels. For shallow water, the flows of liquids on solids are mainly represented by using the Navier-Stokes equations associated with adherence conditions at the walls. Recent experiments in nanofluidics seem to prove, also for liquids, that at nanoscales corresponding to sap layers at very high tree levels, the conditions of adherence are disqualified [START_REF] De Gennes | On the fluid/wall slippage[END_REF][START_REF] Tabeling | Microfluidics[END_REF]. With the aim of explaining experimental results, we reconsider the fluids as media whose motions generate slips along the walls; so, we can draw consequences differing from results of classically adopted models as reconsidered in [START_REF] Gouin | Dynamics of liquid nanofilms[END_REF]. The new model we are presenting reveals an essential difference between the flows of microfluidics and those of nanofluidics. In the latter, simple laws of scales cannot be only taken anymore into account. The transpiration in the leaves induces a variation of the sap layer thickness in microtubes. Consequently, the gradient of thickness along microtubes creates a gradient of disjoining pressure which induces driving forces along the layer. For thin layers, the sap flow depending on the variations of the layer thickness can be adapted to each level of leaves following the tree requirement. This is an important understanding why the flow of sap can be non negligible at a level corresponding to the top of the tallest trees. Moreover, we notice that the stability criterium of the flow issued from the equation of motion fits with the results of Derjaguin's school [START_REF] Derjaguin | Surfaces Forces[END_REF].

A study of inhomogeneous fluids near a solid wall

In this section, we recall the main results presented in [START_REF] Gouin | A new approach for the limit to tree height using a liquid nanolayer model[END_REF][START_REF] Gouin | Liquid nanofilms. A mechanical model for the disjoining pressure[END_REF]. Thanks to these results, in Section 3, we shall consider sap layers of the highest trees with a thickness of some nanometers only. The density-functional of an inhomogeneous fluid in a very thin isothermal layer domain (O) of wall boundary (S) and liquid-vapor interface (Σ) was chosen in the form:

F = (O) ε dv + (S) φ ds + (Σ) ψ ds. (1)
• The first integral is associated with a square-gradient approximation when we introduce a specific free energy of the fluid at a given temperature T as a function of density ρ and β = (grad ρ) 2 such as [START_REF] Widom | What do we know that van der Waals did not know?[END_REF][START_REF] Rowlinson | Molecular Theory of Capillarity[END_REF]:

ρ ε = ρ α(ρ) + λ 2 (grad ρ) 2 ,
where term (λ/2) (grad ρ) 2 is added to the volume free energy ρ α(ρ) of a compressible fluid and scalar λ is assumed to be constant at a given temperature [START_REF] Rocard | Thermodynamique[END_REF]. Specific free energy α enables liquid and vapor bulks to be continuously connected and the pressure P (ρ) = ρ 2 α ′ ρ (ρ) is similar to van der Waals one.

• For a plane solid wall (S), the solid-liquid surface free energy is in the form [START_REF] Gouin | Energy of interaction between solid surfaces and liquids[END_REF][START_REF] Weiss | Theoretical description of the adsorption and the wetting behavior of alkanes on water[END_REF]:

φ(ρ) = -γ 1 ρ + 1 2 γ 2 ρ 2 . ( 2
)
Here ρ denotes the fluid density value at surface (S); constants γ 1 , γ 2 and λ are generally positive and given by the mean field approximation:

γ 1 = πc ls 12δ 2 m l m s ρ sol , γ 2 = πc ll 12δ 2 m 2 l , λ = 2πc ll 3σ l m 2 l , (3)
where c ll and c ls are two positive constants associated with Hamaker constants; σ l and σ s denoting fluid and solid molecular diameters, δ = 1 2 (σ l + σ s ); m l , m s denote masses of fluid and solid molecules; ρ sol is the solid density.

• For the plane liquid-vapor interface (Σ) the surface free energy ψ is reduced to [START_REF] Gouin | A new approach for the limit to tree height using a liquid nanolayer model[END_REF][START_REF] Gouin | Liquid nanofilms. A mechanical model for the disjoining pressure[END_REF]:

ψ(ρ) = γ 4 2 ρ 2 ,
where ρ is the density of the liquid bounding the interface and γ 4 is associated with the interfacial thickness of the order of the fluid molecular diameter (γ 4 ≃ γ 2 ).

In case of equilibrium, functional ( 1) is stationary and yields the equation of equilibrium and boundary conditions [START_REF] Gouin | Boundary conditions for a capillary fluid in contact with a wall[END_REF][START_REF] Pismen | Disjoining potential and spreading of thin liquid layers in the diffuse-interface model coupled to hydrodynamics[END_REF].

Equation of equilibrium

The equation of equilibrium is [START_REF] Gouin | Boundary conditions for a capillary fluid in contact with a wall[END_REF][START_REF] Gouin | Utilization of the second gradient theory in continuum mechanics to study the motion and thermodynamics of liquid-vapor interfaces[END_REF]:

div σ + ρ g i = 0 , (4)
where σ = -ρ 2 ε ′ ρ -ρ div (λ grad ρ) 1 -λ grad ρ ⊗ grad ρ, g is the acceleration of gravity and i, of coordinate x, is the downward direction. Let us consider an isothermal vertical film of liquid; then in orthogonal system, the coordinate z being external and normal to the flat vertical solid wall, spatial density derivatives are negligible except in direction of z. In the complete liquid-vapor layer (we call interlayer ) and along direction z, Eq. ( 4) yields a constant value at each level x:

P (ρ) + λ 2 dρ dz 2 -λ ρ d 2 ρ dz 2 = P v bx ,
where P v bx = P (ρ v bx ) denotes the pressure in the vapor bulk of density ρ v bx bounding the liquid layer at level x. In the fluid, Eq. ( 4) can be written [START_REF] Gouin | Utilization of the second gradient theory in continuum mechanics to study the motion and thermodynamics of liquid-vapor interfaces[END_REF]:

grad (µ o -λ ∆ρ -g x) = 0, (5)
where µ o is the chemical potential (at temperature T ), chosen null in the liquid and vapor bulks of phase equilibrium densities ρ l and ρ v , respectively; ∆ denotes the Laplacian. Thanks to Eq. ( 5), we obtain in all the fluid and not only in the interlayer :

µ o (ρ) -λ∆ρ -g x = µ o (ρ b ), (6)
where µ o (ρ b ) is the chemical potential value of a liquid mother bulk of density

ρ b such that µ o (ρ b ) = µ o (ρ v b )
, where ρ v b is the density of the vapor bulk bounding the layer at level x = 0 [START_REF] Derjaguin | Surfaces Forces[END_REF]. Equation ( 6) is valid in the interlayer and yields the equation of density profile:

λ d 2 ρ dz 2 = µ bx (ρ), with µ bx (ρ) = µ o (ρ) -µ o (ρ bx ). ( 7
)
where ρ bx is the liquid mother bulk density at level x.

The disjoining pressure for vertical liquid films

The disjoining pressure at level x can be written as [START_REF] Derjaguin | Surfaces Forces[END_REF]:

Π = P v bx -P bx ,
where P bx = P (ρ bx ). At a given temperature T , Π is a function of ρ bx . The reference chemical potential linearized near ρ l is µ o (ρ) = (c 2 l /ρ l )(ρ -ρ l ) where c l is the isothermal sound velocity in liquid bulk ρ l at temperature T [START_REF] Gouin | Bubble number in a caviting flow[END_REF]. In the liquid part of the liquid-vapor film, Eq. ( 7) writes:

λ d 2 ρ dz 2 = c 2 l ρ l (ρ -ρ b ) -g x ≡ c 2 l ρ l (ρ -ρ bx ) with ρ bx = ρ b + ρ l c 2 l g x. (8)
At level x = 0, the liquid mother bulk density is closely equal to ρ l and because of Rel. ( 8), Π can be considered as a function of x [START_REF] Gouin | A new approach for the limit to tree height using a liquid nanolayer model[END_REF]:

Π(x) = -ρ l g x 1 + g x 2 c 2 l . (9)
Now, we consider a film of thickness h x at level x; the density profile in the liquid part of the liquid-vapor film is solution of the system:

         λ d 2 ρ dz 2 = c 2 l ρ l (ρ -ρ bx ), with λ dρ dz | z=0 = -γ 1 + γ 2 ρ | z=0 and λ dρ dz | z=hx = -γ 4 ρ | z=hx .
Let quantities τ , d and γ 3 be defined as:

τ ≡ 1 d = c l λρ l and γ 3 ≡ λτ , (10)
such that d is a reference length. Due to the fact that ρ bx ≃ ρ b ≃ ρ l [START_REF] Derjaguin | Surfaces Forces[END_REF], the disjoining pressure reduces to [START_REF] Gouin | Liquid nanofilms. A mechanical model for the disjoining pressure[END_REF]:

Π(h x ) = 2 c 2 l ρ l (γ 1 -γ 2 ρ l )(γ 3 + γ 4 )e hxτ + (γ 2 -γ 3 )γ 4 ρ l × (γ 2 + γ 3 )γ 4 ρ l -(γ 1 -γ 2 ρ l )(γ 3 -γ 4 )e -hxτ [(γ 2 + γ 3 )(γ 3 + γ 4 )e hxτ + (γ 3 -γ 4 )(γ 2 -γ 3 )e -hxτ ] 2 . ( 11
)
The disjoining pressure of the mixture of liquid and perfect gas is the same than for a single van der Waals fluid and calculations and results are identical to those previously obtained [START_REF] Gouin | A new approach for the limit to tree height using a liquid nanolayer model[END_REF].

Water wetting a vertical plane wall of xylem

Our aim is to point out an example such that previous results provide a value of maximum height for a vertical water film wetting a plane wall of xylem and to estimate the sap layer thickness at this altitude. As proved by Derjaguin et al in [START_REF] Derjaguin | Surfaces Forces[END_REF], (Chapter 2), the Gibbs free energy per unit area G can be expressed as a function of h x :

G(h x ) = +∞ hx Π(h) dh,
where h x = 0 is associated with the dry wall in contact with the vapor bulk and h x = +∞ is associated with a wall in contact with the liquid bulk. The spreading coefficient is S = γ SV -γ SL -γ LV , where γ SV , γ SL , γ LV are the solid-vapor, solid-liquid and liquid-vapor free energies per unit area of interfaces, respectively. The energy of the liquid layer per unit area can be written

E = γ SL + γ LV + G(h x ).
The coexistence of two film segments with different thicknesses is a phenomenon interpreted with the equality of chemical potentials and the equality surface tensions of the two films. A spectacular case corresponds to the coexistence of a liquid film of thickness h p and a dry solid wall associated with h x = 0. The film is the so-called pancake layer corresponding to the condition [START_REF] Derjaguin | Surfaces Forces[END_REF][START_REF] De Gennes | Wetting: statics and dynamics[END_REF]:

G(0) = G(h p ) + h p Π(h p ).
Liquid films of thickness h x > h p are stable and liquid films of thickness h x < h p are metastable or unstable. For a few nanometer range, the film thickness is not exactly h x ; we must add the thickness estimated at 2 σ l of the liquid part of the liquid-vapor interface bordering the liquid layer and the film thickness is e x ≈ h x + 2 σ l [START_REF] Rowlinson | Molecular Theory of Capillarity[END_REF][START_REF] Rocard | Thermodynamique[END_REF]. When h x = 0 (corresponding to the dry wall), the value of G is the spreading coefficient S. Point P associated with the pancake layer is observed on the curve to be closely an inflexion point of graph Π(h x ) [START_REF] Gouin | A new approach for the limit to tree height using a liquid nanolayer model[END_REF]. To obtain the pancake thickness corresponding to the smallest film thickness, we draw the graphs of Π(h x ) and G(h x ) when h x ∈ [(1/2) σ l , ℓ], where ℓ is a distance of few tens of Amgström.

For the numerical calculations, we considered water at T = 20 • Celsius wetting a wall in xylem. The experimental estimates of coefficients are obtained in c.g.s. units [START_REF] Israelachvili | Intermolecular Forces[END_REF][START_REF] Gouin | Liquid nanofilms. A mechanical model for the disjoining pressure[END_REF][START_REF]Handbook of Chemistry and Physics[END_REF]: ρ l = 0.998, c l = 1.478×10 5 , c ll = 1.4×10 -58 , σ l = 2.8×10 -8 (2.8 ångströms), m l = 2.99 × 10 -23 . From Rel. (3), we deduce λ = 1.17 × 10 -5 , γ 2 = γ 4 = 54.2. From Rel. (10), we get γ 3 = 506, d = 2.31 × 10 -8 . We consider a material such that the Young angle between the liquid-vapor interface and the solid material surface is θ ≈ 50 • . This Young angle is an arithmetic average of Young angles for water wetting different xylem walls [START_REF] Mattia | Review: static and dynamic behavior of liquids inside carbon nanotubes[END_REF]. The coefficient c ls is obviously not given in the physical tables and γ 1 cannot be obtain from Rel. (3); due to Rel. (2), we immediately get the unknown coefficient γ 1 in Rel. ( 11): γ 1 ≈ 75.

In Fig. 1 -left graph, we present the disjoining pressure graph Π(h x ). The physical part of disjoining pressure graph corresponding to ∂Π/∂h x < 0 is a plain line and is associated with thickness liquid layer of several molecules. The dashed line has no real existence. In Fig. 1 -right graph, we present the free energy graph G(h x ). Due to h x > (1/2) σ l , it is not possible to obtain the limit point W corresponding to the dry wall. This point is obtained by an interpolation associated with the concave part of the G-curve. Point P follows from the drawing of the tangent line issued from W to the G-curve. The limit of the film thickness is associated to the pancake thickness e p ≈ h p + 2 σ l when the liquid film coexists with the dry wall. The reference length d is of the same order than σ l , σ s and δ and is a good length order for very thin films. The total pancake thickness e p is of one nanometer order corresponding to a good thickness value for a highenergy surface [START_REF] De Gennes | Capillarity and Wetting Phenomena: Drops, Bubbles, Pearls, Waves[END_REF]. We deduce S ≈ 40 cgs. However, crude sap is not pure water. Its liquid-vapor surface tension has a lower value than surface tension of pure water (72 cgs at 20 • C) and it is possible to obtain the same spreading coefficients with less energetic surfaces. When |x| is of some hundred meters, Eq. (9) yields Π(x) ≃ -ρ l g x. The maximum of altitude |x M | corresponds to the pancake layer. To this altitude, we add 20 meters corresponding to the ascent of sap due to capillarity and osmotic pressure and we obtain a film height of 140 meters which is of the same order than the topmost trees. We also note the important result: in the trees, the thickness of the layer is of some nanometers at high level.

Dynamics of liquid nanofilms

The idea of xylem tubes completely filled with sap induces that the flow of liquid along vessels can be compared with a flow through capillaries [START_REF] Zimmermann | Xylem structure and the ascent of sap[END_REF]. The flow rate through capillary tubes is proportional to the applied pressure gradient and the hydraulic conductivity. Parabolic flow causes the flow rate to be proportional to the fourth power of capillary radius [START_REF] Batchelor | An introduction to fluid dynamics[END_REF]. One of the most important design requirements is that vapor blockage does not happen in the stem. Consequently, to be efficient for the transport of sap, the xylem tube radius must be as larger as possible, which is not the case.

The dynamics of such liquid nanofilms is always studied in isothermal case.

Our model of sap nanolayer implies different consequences:

-The classical model by Navier-Stokes is unable to describe fluid motions in very thin films, -The notion of surface free energy of a sharp interface separating gas and liquid layer must be reconsidered, -The equation of fluid motions along the nanofilm is obtained by adding the forces of viscosity to the conservative forces, -The evolution equation of the film thickness takes account of the variation of the disjoining pressure along the layer. At equilibrium, the different fluid quantities are ρ bx , ρ v bx , P bx , P v bx , h x . . . For the motion, the corresponding quantities are denoted by ρ * bx , ρ * v bx , P * bx , P * v bx , h * x . . . When the liquid layer thickness is small with respect to transverse dimensions of the wall, it is possible to simplify the Navier-Stokes equation which governs the flow of classical viscous fluids. When h ≪ L, where L is the wall transversal characteristic size, i) The velocity component along the wall is large with respect to the normal velocity component which can be neglected, ii) The velocity vector mainly varies along the direction orthogonal to the wall and it is possible to neglect velocity derivatives with respect to coordinates along the wall compared to the normal derivative, iii) The pressure is constant in the direction normal to the wall. It is possible to neglect the inertial term when Re ≪ L/h, where Re is the Reynolds number of the flow. Equation of hydrodynamics is not valid in a liquid nanolayer because the fluid is heterogeneous and the liquid stress tensor is not anymore scalar. However, it is possible to adapt the results obtained for viscous flows to motions in liquid nanolayers: due to ǫ = h/L ≪ 1, we are also in the case of long wave approximation. We denote the velocity by V = (u, v, w) where (u, v) are the tangential components. Due to the fact that e = sup (|w/u| , |w/v|) ≪ 1, we also are in the case of approximation of lubrication. The main parts of terms associated with second derivatives of liquid velocity components correspond to ∂ 2 u/∂z 2 and ∂ 2 v/∂z 2 . The density is constant along each stream line ( ρ = 0 ⇐⇒ div V = 0) and isodensity surfaces contain the trajectories. Then, ∂u/∂x, ∂v/∂y and ∂w/∂z have the same order of magnitude and ǫ ∼ e. As in [START_REF] Rocard | Thermodynamique[END_REF], we assume that the kinematic viscosity coefficient ν = κ/ρ only depends on the temperature. In motion equation, the viscosity term is

(1/ρ) div σ v = 2ν [ div D + D grad { Ln (2 κ)} ] ,
where D is the velocity deformation tensor and D grad{Ln (2 κ)} is negligible with respect to div D . In both lubrication and long wave approximations, the liquid nanolayer motion verifies [START_REF] Gouin | Utilization of the second gradient theory in continuum mechanics to study the motion and thermodynamics of liquid-vapor interfaces[END_REF]:

Γ + grad[ µ o (ρ) -λ ∆ρ ] = ν ∆V + g i with ∆V ≃ ∂ 2 u ∂z 2 , ∂ 2 v ∂z 2 , 0 .
This equation corresponds to the equation of equilibrium (4) with addition of inertial (acceleration) term Γ and viscous term ν ∆V. In approximation of lubrication, the inertial term is neglected:

grad[ µ o (ρ) -λ ∆ρ ] = ν ∆V + g i. (12)
Equation ( 12) separates into tangential and normal components to the solid wall.

-The normal component -following z -of Eq. ( 12) writes in the same form than for equilibrium:

∂ ∂z [ µ o (ρ) -λ ∆ρ ] = 0,
and consequently,

µ o (ρ) -λ ∆ρ = µ o (ρ * bx ), (13)
where ρ * bx is the dynamical liquid mother bulk density at level x (different from the liquid bulk density ρ l of the plane interface at equilibrium and also from ρ bx , liquid mother bulk density at level x at equilibrium). To each density ρ * bx is associated a liquid nanolayer thickness h * x . We can write

µ o (ρ b ) -λ ∆ρ = η(h * x ), where η is such that η(h * x ) = µ o (ρ * bx ).
-For one-dimensional motions colinear to the solid wall (direction i and velocity u i), by taking account of Eq. ( 13), the tangential component of Eq. ( 12) yields:

i . grad µ o (ρ * bx ) = ν ∂ 2 u ∂z 2 + g ,
which is equivalent to:

∂µ o (ρ * bx ) ∂ρ * bx ∂ρ * bx ∂x = ν ∂ 2 u ∂z 2 + g. ( 14
)
For the most practical situations, simple fluids slip on a solid wall only at a molecular level [START_REF] Chuarev | Thin liquid layers[END_REF] and consequently, the kinematic condition at solid walls is the adherence condition (z = 0 ⇒ u = 0). Nevertheless, with water flowing on thin hydrophobic nanolayers, there are some qualitative observations for slippage. With water flowing in thin, hydrophobic capillaries, there are also some early qualitative evidences for slippage [START_REF] Blake | Slip between a liquid and a solid -D.M. Tolstoi (1952) theory reconsidered[END_REF]. De Gennes said: "the results are unexpected and stimulating and led us to think about unusual processes which could take place near a wall. They are connected with the thickness h of the film when h of an order of the mean free path" [START_REF] De Gennes | On the fluid/wall slippage[END_REF].

Recent papers in nonequilibrium molecular dynamics simulations of three dimensional micro-Poiseuille flows in Knudsen regime reconsider microchannels: the influence of gravity force, surface roughness, surface wetting condition and wall density are investigated. The results point out that the no-slip condition can be observed for Knudsen flow when the surface is rough. The roughness is a dominant parameter as far as the slip of fluid is concerned. The surface wetting condition substantially influences the velocity profiles [START_REF] Tabeling | Microfluidics[END_REF].

In fluid/wall slippage, the condition at solid wall writes:

u = L s ∂u ∂z at z = 0,
where L s is the so-called Navier length. The Navier length is expected to be independent of h and may be as large as a few microns [START_REF] Tabeling | Microfluidics[END_REF]. At the liquid-vapor interface, we also assume that vapor viscosity stress is negligible; from the continuity of fluid tangential stress through a liquid-vapor interface, we get:

∂u ∂z = 0 at z = h * x .
Consequently, Eq. ( 14) implies:

ν u = ∂µ o (ρ * bx ) ∂ρ * bx ∂ρ * bx ∂x -g 1 2 z 2 -h * x z -L s h * x .
The mean spatial velocity u of the liquid in the nanolayer at level x is:

u = 1 h * x h * x o u dz.
Previous computations yield:

ν u = -h * x h * x 3 + L s grad µ o (ρ * bx ) -g i with u = u i .
Let us note that:

∂µ o (ρ * bx ) ∂x = ∂µ o (ρ * bx ) ∂ρ * bx ∂ρ * bx ∂h * x ∂h * x ∂x ≡ 1 ρ * bx ∂P (ρ * bx ) ∂ρ * x ∂ρ * bx ∂h * x ∂h * x ∂x .
Due to the fact the vapor bulk pressure P * v bx is constant along the xylem tube, by using relation Π(h * x ) = P * v bx -P * bx , we get along the flow motion:

∂µ o (ρ * bx ) ∂x = - 1 ρ * bx ∂Π(h * x ) ∂h * x ∂h * x ∂x and consequently, χ * bx u = h * x h * x 3 + L s [ grad Π(h * x ) + g i] , (15)
where χ * bx = ρ * bx ν is the liquid kinetic viscosity. Consequently, the mean liquid velocity is driven by the variation of the disjoining pressure along the solid wall and the film thickness. Equation ( 15) differs from the classical film hydrodynamic one. Indeed, for a classical thin liquid film, the Darcy law is u = -K(h) grad p, where p is the liquid pressure and K(h) is the permeability coefficient. In Eq. ( 15), the sign is opposite and the liquid pressure is replaced by the disjoining pressure. We note that χ * bx ≃ χ , where χ is the liquid kinetic viscosity in the liquid bulk at phase equilibrium. Moreover, when h * x /L s ≪ 1,

χ u = h * x L s [ grad Π(h * x ) + g i] ,
which is strongly different from the case L s = 0 corresponding to the adherence condition:

χ u = h * 2 x 3 [ grad Π(h * x ) + g i] .
The mass equation averaged over the liquid depth is:

∂ ∂t h * x 0 ρ dz + div h * x 0 ρ u dz = 0.
Since the variation of density is small in the liquid nanolayer, the equation for the free surface is:

dh * x dt + h * x div u = 0. (16)
By replacing ( 15) into ( 16) we finally get:

∂h * x ∂t + 1 χ div h * 2 x h * x 3 + L s grad Π(h * x ) + g i = 0. ( 17
)
It is easy to verify that Eq. ( 17) is a non-linear parabolic equation.

If ∂Π(h * x )/∂h * x < 0 the flow is stable. This result is in accordance with the static criterium of stability for thin liquid layers [START_REF] Derjaguin | Surfaces Forces[END_REF]. When L s = 0 we notice the flow is multiplied by the factor 1 + 3L s /h * x . For example, when h * x = 3 nm and L s = 100 nm which is a Navier length of small magnitude with respect with experiments, the multiplier factor is 10 2 ; when L s is 7 µm as considered in [START_REF] Tabeling | Microfluidics[END_REF], the multiplier factor is 10 4 , which seems possible in nanotube observations [START_REF] Mattia | Review: static and dynamic behavior of liquids inside carbon nanotubes[END_REF]. Equation ( 15) yields the flow rate per unit of length of xylem tubes. We may remark that Eq. ( 15) is mainly valid at the top of highest trees where the xylem tube network is strongly ramified. A main difference between Poiseuille flow and motion in a thin film is the versatility of the liquid layer flow with respect to Poiseuille'. An hydraulic Poiseuille flow is very rigid due to the liquid incompressibility, the pressure effects are fully propagated in all the tube. For a thin layer flow, the flow rate can increase or decrease due to the spatial derivative of h *

x and strongly depends on the locally defined disjoining pressure. Trees can adapt the disjoining pressure effects by opening or closing the stomatic cells with the object of changing the evaporation in its leaves so that the bulk pressure in micropores can be negative and the transport of water can be differently dispatched in the tree parts; this seems an important aspect of the model.

Discussion and conclusions

Our model is essentially different from the cohesion-tension theory; it allows a new explanation of biofluidics by using methods of the Russian school of Derjaguin and using non-adherence conditions for nanofluid flows at xylem walls. It precisely gives the possible height of the tallest trees. In the very thin layers, we obtain rates of flows strongly larger than those obtained with traditional Navier-Stokes models. The explanation of irrigation of the leaves in the tallest trees supplementary justifies our model of biomimetism. The motor of the sap motion is induced by the transpiration across micropores located in tree leaves [START_REF] Zimmermann | Xylem structure and the ascent of sap[END_REF]. It is natural to conjecture that the diameters of xylem tubes must be the result of a competition between evaporation in tubes which reduces the flow of sap and the flux of transpiration in micropores inducing the motion strength. We notice that the negative pressure only appears in the liquid mother bulk. The microlimbs, micropores and stomates have a diameter a little smaller than the bubble size considered in the introduction; consequently, they can be filled without any cavitation by the liquid mother bulk at a suitable negative pressure associated with the height of the tallest trees. It is interesting to note that if we switch the microtube surfaces to wedge geometry or to corrugated surface, it is much easier to obtain the complete wetting requirement; thus, plants can avoid having very high energy surfaces. Nonetheless, they are still internally wet if crude sap flows pass through wedge shaped corrugated pores; this fact also answers to questions in [START_REF] Zimmermann | Xylem structure and the ascent of sap[END_REF]. The wedge does not have to be perfect on the nanometer scale to significantly enhance the amount of liquid flowing at modest pressures corresponding to nanosized planar films. It is bound to improve on the calculation because it enhances the surface to volume ratio. In such a case, we remark that the wall boundary can always be considered as a plane surface endowed with an average surface energy as in Wenzel's formula [START_REF] Wenzel | Surface roughness and contact angle[END_REF]. Finally, it will be interesting to confirm our theoretical predictions with additional experimental data.

Fig. 1 .

 1 Fig. 1. Left graph: Π(h x )-graph. The unit of x-axis is d = 2.31 × 10 -8 cm; the unit of y-axis is one atmosphere. Right graph: G(h x )-graph. The unit of x-axis is d = 2.31 × 10 -8 cm ; the unit of y-axis is one cgs unit of surface tension.