Tidally driven pore water exchange within offshore intertidal sandbanks: Part II numerical simulations

Estuarine, Coastal and Shelf Science 80 (2008) 472–482

Contents lists available at ScienceDirect

Estuarine, Coastal and Shelf Science journal homepage: www.elsevier.com/locate/ecss

Tidally driven pore water exchange within offshore intertidal sandbanks: Part II numerical simulations B. Gibbes a, *, C. Robinson a, L. Li a, b, D. Lockington a, H. Li c a

Environmental Engineering Division, School of Engineering, The University of Queensland, St. Lucia QLD 4072, Australia Centre for Eco-Environmental Modeling, State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098, P.R.China c School of Environmental Studies, China University of Geosciences, Wuhan 430074, P.R.China b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 20 July 2007 Accepted 28 August 2008 Available online 25 September 2008

Field measurements presented by [Gibbes, B., Robinson, C., Li, L., Lockington, D.A., Carey, H., 2008. Tidally driven pore water exchange within offshore intertidal sandbanks: Part I Field measurements. Estuarine, Coastal and Shelf Science 79, pp. 121–132.] revealed a tidally driven pore water flow system within an offshore intertidal sandbank in Moreton Bay, Australia. The field data suggested that this flow system might be capable of delivering nutrients, and in particular bio-available iron, across the sediment–water interface. Bio-available iron has been implicated as a key nutrient in the growth of the toxic marine cyanobacteria Lyngbya majuscula and therefore this pore water exchange process is of interest at sites where L. majuscula blooms have been observed. In this study two-dimensional numerical simulations were used in conjunction with hydraulic data from field measurements to further investigate the tidally induced pore water flow patterns. Simulation results generally showed good agreement with the field data and revealed a more complex residual pore water flow system in the sandbank than shown by the field data. The flow system, strongly influenced by the geometry of the sandbank, was characterized by two circulation cells which resulted in pore water discharge at the bank edge and also to a permanently ponded area within the sandbank interior. Simulated discharge volumes in these two zones were in the order of 0.813 m3 and 0.143 m3 per meter width (along shore) of sandbank per tidal cycle at the bank edge and sandbank interior respectively. Transit times of pore water circulating through these cells were found to range from z 17 days to > 60 years with an average time of 780 days. The results suggest that the tidally driven flow systems might provide a mechanism for transport of bio-available iron across the sediment–water interface. This flow could constitute a previously unrecognized source of bio-available iron for L. majuscula blooms in the Bay. Ó 2008 Elsevier Ltd. All rights reserved.

Keywords: intertidal system pore water exchange numerical model offshore sandbank Lyngbya majuscula Moreton Bay

1. Introduction The action of tides can induce significant pore water exchange between intertidal sediments and the adjacent coastal water body in nearshore salt marsh environments (Ursino et al., 2004; Gardner, 2005; Li et al., 2005; Billerbeck et al., 2006; Gardner and Wilson, 2006; Wilson and Gardner, 2006). In these environments water is drained through the edge of the salt marsh during the ebb tide and is then replaced by water infiltrating through the upper salt marsh surface as the marsh is overtopped (submerged) on the flood tide. Recent field data has shown that a similar tidal exchange mechanism also operates in offshore intertidal sandbanks (Gibbes et al., 2008). In these sandbanks the tidal exchange mechanism is likely to be capable of delivering reduced chemical species across the

* Corresponding author. E-mail address: [email protected] (B. Gibbes). 0272-7714/$ – see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ecss.2008.08.021

sediment–water interface, which might constitute a potentially significant input of nutrients to coastal water bodies. This potential nutrient delivery mechanism is of particular importance for offshore intertidal sandbanks in Moreton Bay, Australia, where blooms of the toxic marine cyanobacterium Lyngbya majuscula (family Oscillatoriaceae, referred to herein as L. majuscula) are observed. Bio-available forms of iron (Fe(II)) have been implicated as a key nutrient in L. majuscula growth (Albert et al., 2005; O’Neil and Dennison, 2005; Watkinson et al., 2005; Ahern et al., 2006, 2007). While the data presented by Gibbes et al. (2008) provides new information on the tidally induced flow system in the sandbank, limited placement of pressure transducers constrained the development of a more complete picture of pore water exchange between the sandbank and the Bay. As with nearshore salt marshes, numerical models of pore water flow can be used to further investigate this tidally induced flow system. Development of models that can effectively simulate such a flow system is challenging. As highlighted by Gardner and Wilson

B. Gibbes et al. / Estuarine, Coastal and Shelf Science 80 (2008) 472–482

(2006) pore water flows in nearshore salt marshes are driven by complex interactions among semi-diurnal tidal forcing, diurnal evapotranspiration, episodic rain events, saturated and unsaturated flow as well as the development of seepage faces. Offshore intertidal sandbanks are in many ways similar to nearshore salt marsh systems. They are hydrologically similar to the ‘‘marsh islands’’ described by Novakowski et al. (2004) in that there is no terrestrially derived or ‘‘fresh’’ groundwater flow and hence consideration of density-driven flows due to the interaction of seawater and low-salinity groundwater is not required. In these situations tidal exchange can become the dominant pore water exchange process. Recent model development for nearshore salt marsh systems (Harvey et al., 1987; Ursino et al., 2004; Gardner, 2005; Li et al., 2005; Gardner and Wilson, 2006) provides a useful guide to the techniques that can be applied to offshore intertidal sandbanks. To the authors’ knowledge, modelling of pore water flow in offshore intertidal sandbank systems has not been reported in the literature. The common occurrence of these intertidal sandbanks in marine tidal delta and some fluvial delta zones of many estuarine environments (Roy et al., 2001) combined with their large areal extent suggests that they are likely to play an important role in the functioning and structure of coastal ecosystems. Furthermore much of the existing model development for nearshore salt marsh systems has, for various reasons, focussed on idealized systems (rather than a specific environmental system or field site). Reported comparisons between model predictions and field data for a specific site are rare (e.g., Harvey et al. (1987)). In this study numerical simulations are conducted to further investigate the tidally induced pore water flow patterns at the field site described by Gibbes et al. (2008). Through comparison with collected field data, the simulations are used to provide further insights into the pore water flow system in the offshore sandbanks. Particular emphasis is placed on the rates of pore water discharge across the sediment–water interface. It is recognized that coupling of the flow model to a reactive transport model would allow further investigation of the influence of the tidal flow systems on the fate and fluxes of reactive chemical species, such as bio-available iron. However, incorporation of proper reactive transport processes into the flow model is a challenging task and is thought to be pre-mature at this stage with limited geochemical data. Therefore the focus of this investigation is the pore water flow system in the offshore intertidal sandbank. The findings will also have implications for studies of flows in nearshore salt marshes. 2. Numerical model A model based on the time-dependent pressure form of the Richards equation for variably saturated porous medium flow was used to simulate tidally driven exchange in a two-dimensional, shore normal, vertical section of the offshore sandbank. The governing equation for the model is given by:

" #   ks  C vp V$ kr V p þ rf gD ¼ þ Se S m rf g vt

i h S ¼ Xp ð1  qs Þ þ Xf qs

473

(2)

where qs is the liquid volume fraction [-], Xp and Xf are the compressibility [M1 L T2] of the solid and the fluid respectively. This form of the Richards equation does not incorporate the total stress term discussed by Reeves et al. (2000). The lag and damping effects associated with the storage term under tidal loading were effectively reduced to a negligible level by altering the storage term to minimize the apparent aquifer compression. The alteration of the storage terms in effect assumes that neither the water nor aquifer matrix is compressible and only allows for changes in water content (i.e., degree of saturation, the overwhelmingly dominant contributor to storage). Simple models based on the governing equation above were developed, using the scenario described by Reeves et al. (2000), and confirmed that the influence of compression was negligible for the sand sediment type used in this study. As such it is assumed to be an acceptable representation of the sand sediment type used in this investigation rather than resorting to a fully coupled model for flow and poroelastic deformation of the porous medium. 2.1. Model domain and boundary conditions Data from the sand level surveys at the Moreton Bay field site (refer to Gibbes et al. (2008)) were used to develop a two-dimensional cross-sectional model of the site (Fig. 1). Sand level data was used to establish the geometry of the upper sandbank surface and sandbank edge (M–N). The boundary N–O was added to the model geometry to ensure that the vertical boundary (O–P) was distant from the active areas of the model thus minimizing the effects of this boundary on the flow field. The upper sandbank surfaces (M–N and N–O) were represented by a mixed form, head-dependent flux boundary condition (Cauchy-type condition) that simulated the inflow and outflow across these boundaries caused by the tideinduced pressure changes:

 i  h !k n Vp þ rf gVD ¼ Rb ðpb  pÞ þ rf gðDb  DÞ

m

(3)

! where n is the vector normal to the boundary, k ¼ ks  kr, Rb is the conductance [M1L 2T] of materials between the source and the model domain, pb and Db are the pressure [M L1 T2] and elevation [L] of a distant fluid source. The Cauchy-type boundary condition incorporated a high conductance for the material between the source and the model domain to effectively allow instant transmission of flow through this material. The elevation of the distant fluid source Db was equal to the boundary elevation. In this

(1)

where V is the gradient operator [L1], ks is the intrinsic permeability [L2], m is the fluid viscosity [M L1 T1], kr is the relative permeability [-], p is the fluid pressure [M L1 T2], rf is the fluid density [M L3], g is the magnitude of gravitational acceleration [L T2], D is the vertical elevation [L], C is the specific moisture capacity [L1], Se is the effective saturation [-], S is the storage coefficient [M1L T2] and t is time [T]. The soil moisture characteristics (kr, C and Se) were given by the van Genuchten (1980) retention and permeability relationships. The storage coefficient (S) was defined by:

Fig. 1. Conceptual illustration of the two-dimensional model domain used to simulate a cross-section at the Moreton Bay field site.

474

B. Gibbes et al. / Estuarine, Coastal and Shelf Science 80 (2008) 472–482

application the fluid source is akin to a tidal reservoir located on the boundary. This methodology provided a good approximation of the condition given by the tide. The semi-diurnal tidal signal at the field site was approximated by a sinusoidal function:

h  pi hðxb ; zb ; tÞ ¼ zMSL þ Acos ut  2

(4)

where h(xb, zb, t) is the local tidal head [L] on the boundary at a given time (xb and zb are the coordinates of the boundary), zMSL is the elevation of the mean tide level [L], A is the tidal amplitude [L] and u is the angular frequency [T1] of the tide. The angular frequency of the tide was defined by:

u ¼

2p P

(5)

where P is the tidal period [T]. Table 1 summarises the tidal parameters used in the simulations. The tidal parameters were chosen to approximate the measured tide data at the field site. The relationship between the tidal head h (equation (4)) and the boundary pressure term pb used in equation (3) was given by equation (6).

pb ¼ rf gh

(6)

A Boolean condition was incorporated into the conductance term (Rb) in the Cauchy-type boundary condition (equation (3)) to set the term, and therefore the boundary flux, to zero when the tidal elevation was below the elevation of the boundary (i.e., Rb ¼ 0 when h < zb). That is, the right-hand-side of equation (3) is set to zero when the elevation of the tide is below the boundary. As a result these tidal boundaries (M–N and N–O) permitted flow across them when they were submerged but acted as no-flow boundaries during exposure. This no-flow boundary condition effectively prevented moisture loss from evaporation from being incorporated in the model. It also prevented the simulation of a seepage face that might form on the sandbank edge. This choice of boundary condition was dictated by the software package which was unable to use a mixed or switching boundary condition (Dirichlet-type to Neumann-type) that is commonly employed in similar tidal models (e.g., Gardner, 2005; Gardner and Wilson, 2006; Wilson and Gardner, 2006). It is acknowledged that such a mixed boundary or switching boundary condition can allow for more physically realistic representation of the system. Despite this limitation the software package was still used due to its ability to easily allow irregular topography to be included in the model (a feature that is thought to have a greater influence over pore water flows in this case). Boundary L–P was located at the assumed depth of the relatively impermeable bedrock material and was represented by a no-flow (Neumann-type) boundary condition. The left and right hand side boundaries (M–L and O–P) were assumed to be flow divides and thus were also represented by no-flow (Neumanntype) boundary conditions. It should be noted that the Cauchy-type boundary condition was used instead of a Dirichlet-type (or specified head) boundary condition on the tidal boundaries as the latter condition was found to cause ‘‘artificial wetting’’ along these boundaries. This artificial wetting is a modelling artefact, which results from the pressure

dependence of the van Genuchten (1980) effective saturation term and hence the water saturation. This pressure dependence can result in mass conservation errors in certain situations where there are dynamic boundary conditions (e.g., a tidally influenced system). Where the artefact becomes most noticeable in the current model is at the time immediately before overtopping of the upper sandbank surface on the rising tide. A Dirichlet-type boundary condition will result in artificial wetting as the negative pressure on the boundary approaches zero before overtopping occurs and the effective saturation starts to increase. This increase in effective saturation is independent of the flow system and can lead to ‘‘artificial wetting’’ of the porous medium close to the boundary. This artefact also prevents a two-zone modelling approach (Winter, 1976; Anderson et al., 2002), which is commonly utilized in coastal groundwater applications (Robinson et al., 2007), from being applied to variably saturated models of tidal systems. 2.1.1. Porous medium properties The porous medium in the model domain was assumed to be homogenous and isotropic. The available data on the hydraulic properties of the site suggest that this is a reasonable assumption (Gibbes et al., 2008). Table 2 provides a summary of the sediment hydraulic properties that were used in the numerical simulations. The average hydraulic conductivity (Ks ¼ 2.025  104 m s1) value derived from the bail tests at the field site (Gibbes et al., 2008) was adopted for the numerical model. The hydraulic conductivity (Ks) is related to the intrinsic permeability (ks) used in equation (1) by equation (7).

ks ¼

Ks m rf g

(7)

Additional model runs were conducted using Ks values that were two and four times as large as the field test value. The purpose of performing simulations with different Ks values is discussed in Section 3.1.1. The moisture retention properties of the sediment at the field site, as used in the numerical model, were assumed to be the same as those presented by Carsel and Parrish (1988) for a sand sediment type. 2.2. Numerical scheme Equation (1) was discretized and solved using the FEMLAB 3.1i Multiphysics Finite Element Analysis Software Package (COMSOL, 2004). Free time stepping was implemented to ensure small time steps in order to avoid the global mass balance errors identified by Celia et al. (1990), which are often associated with the pressurebased form of the Richards equation. Simulations were initiated with uniform pressure heads given by a mean sea level of 0 m Australian Height Datum (AHD). This resulted in an initial condition where the sandbank was completely submerged with a completely saturated porous medium. The numerical model was run for several tidal cycles to ensure that the simulation had reached a quasisteady-state (i.e., periodic with respect to time). Checks were performed on all models to verify that the simulation results were independent of the density of the grid used to define the finite element mesh.

Table 2 Soil moisture retention parameters used in numerical simulations

Table 1 Tidal parameters used in numerical simulations Tide type

Mean tide level (ZMSL) [m]

Amplitude (A) [m]

Period (P) [s]

Frequency (u) [Rad s1]

Semidiurnal

0.0

0.774

4.32  104 (12 h)

1.4544  104

Soil type Sand a

Ks [m s1] 2.025  10

4a

qs [m3 m3]

qr [m3 m3]

a [m1]b

n [-]b

l [-]

0.43

0.045

14.5

2.68

0.5

Two additional simulations were run with Ks ¼ 4.050  104 and 8.101 104 m s1 respectively (other parameters values unchanged). b van Genuchten (1980) coefficients.

B. Gibbes et al. / Estuarine, Coastal and Shelf Science 80 (2008) 472–482

3. Results 3.1. Comparison of simulated and field head gradients The performance of the numerical model was assessed by comparing the simulated vertical head gradients with those observed at the field site. Due to the complete saturation of the sandbank at the end of the overtopping period, the model was found to essentially ‘‘reset’’ on every tidal cycle. Therefore only data for a single representative tidal cycle is presented. The following discussion uses the same reference system presented by Gibbes et al. (2008) to describe the location of the monitoring stations. The same convention for describing the head gradient is also used (i.e., a negative gradient indicates upward flow and a positive gradient indicates downward flow). 3.1.1. Bank platform – station C For monitoring station C there was an overall agreement between the simulation results and field data (Fig. 2). Flow at this location was predominantly downward as shown by both the field data and simulation data. The downward flow reflects the drainage characteristics of the sandbank during exposure. The simulated and measured magnitudes of the positive gradient (downward flow) were also generally consistent. The primary difference between the field and simulation results was that the ponding condition and associated positive peak in the vertical head gradient (Fig. 2(a) and (c) at elapsed time z 129.5 h) observed in the field was not replicated by the model. Comparison of the field and simulation data also showed that the model was able to effectively simulate the rapid overtopping of the bank platform (Fig. 2 at elapsed time z 134.5 h). Prior to overtopping of the bank platform negative vertical gradients (indicating upward flow) were evident in both the field and simulation data. This behaviour was consistent with a rising water table due to the tide. As the bank platform was overtopped, water started to infiltrate the sandbank and the gradient changed direction rapidly (to a positive gradient). The replication of the overtopping event by the model was thought to be particularly

475

important as the rapid overtopping of the bank platform on the rising tide is the key driver for asymmetric flows through the sandbank and hence the tidally driven exchange. The results also indicated that the model based on the Ks value derived from field tests under-predicted the lowering of the water table during the exposure period by approximately 0.3 m (Fig. 2(a) cf. 2(b)). To investigate the influence of Ks values two additional simulations were performed with the hydraulic conductivity value multiplied by a factor of two and four (i.e., Ks ¼ 4.05  104 m s1 and Ks ¼ 8.101 104 m s1 respectively). The results of these additional simulations revealed that with higher Ks values the model predicted a more pronounced water table lowering but not to the same extent as observed in the field data (approximately 0.1 m under-prediction for Ks ¼ 8.101 104 m s1). The different Ks values did not generally have an effect on the vertical head gradients. Simulations that incorporated different Ks values produced similar vertical head gradient patterns. Some variation in the magnitude of the positive (downward flow) vertical head gradient was observed during the overtopping event. The positive vertical head gradient increased with increasing Ks values. An increase by approximately a factor of 2.5 was observed between the simulation which used the average field Ks values (Ks ¼ 2.025  104 m s1) and that using Ks ¼ 8.101 104 m s1 (maximum positive vertical head gradients of 0.022 and 0.057 respectively). The maximum vertical head gradient measured at station C at the field site during this overtopping event (0.022) was the same as that simulated by the model which used the average field Ks values. 3.1.2. Middle intertidal bank edge – stations F and G At station F (Fig. 3 (a–d)), there was good agreement between the field data and model predictions during the early and middle stages of sandbank exposure. As the falling shoreline crossed station F on the ebbing tide (Fig. 3(a–d), elapsed time z 129 h) both the model and the field data showed a distinct negative gradient (upward flow) that was quickly followed by a positive gradient (downward flow). Note that the shoreline is defined herein as the intersection between the sea level and sandbank edge. The

Fig. 2. Data for monitoring station C (Middle bank platform) comparing (a) the observed head and (c) vertical head gradient data with the (b) simulated head and (d) simulated vertical head gradient data. The solid black line (–) in (a) and (b) represents the tidal water level while the short-spaced dashed line (- - -) is the head within the sandbank. The horizontal dashed line (– –) represents the sandbank elevation at the monitoring station. Negative gradient values indicate upward flow.

476

B. Gibbes et al. / Estuarine, Coastal and Shelf Science 80 (2008) 472–482

magnitudes of the negative vertical head gradient given by the field data and simulation data were similar (z 0.03, Fig. 3(c) and (d)). During the middle stages of bank exposure both the field and simulation data showed a small positive gradient (downward flow) as the local flow over this period was expected to be predominantly horizontal. Analysis of the horizontal gradient data from the model confirmed such a horizontal flow, with the same order of magnitude as that of the peak vertical flow (z0.03–0.04). The simulated head gradient behaved differently from the field measurement during the overtopping event (Fig. 3(c) and (d) elapsed time z133.5 h). While only positive gradients (downward flow) were evident in the field data during overtopping, the model predicts a pulse of alternating gradient similar to that associated with the passage of the shoreline across the monitoring point on

the ebbing tide (discussed above). The cause of this difference was not readily apparent. Comparison of the gradients from models with different Ks values showed the models with larger Ks values replicated more closely the later stage of the overtopping event where the gradient is positive (downward flow). The large Ks values did not enable the models to replicate the field data during the initial stages of the overtopping event. At monitoring station G the simulation results agreed with the field data around the early stage of exposure (elapsed time z129.5 h, Fig. 3(e–h)). This period was characterized by a strong negative pulse in the head gradient data (upward flow). The peak gradient coincided with the shoreline crossing the monitoring location. The magnitude of the vertical head gradients produced by the numerical simulation was consistent with that observed at the

Fig. 3. Data for monitoring stations F and G (middle intertidal bank edge) comparing the observed head ((a) and (e)) and vertical head gradient ((c) and (g)) with the simulated head ((b) and (f)) and simulated vertical head gradient ((d) and (h)). All other presentation conventions are the same as in Fig. 2.

B. Gibbes et al. / Estuarine, Coastal and Shelf Science 80 (2008) 472–482

field site. Once the falling shoreline passed the monitoring point the field data exhibited a positive vertical gradient (downward flow). The vertical flow then became negligible (zero gradient) for most of the exposure period until the subsequent overtopping. The simulation data did not replicate this positive gradient. The simulated gradient instead diminished following the downward peak gradient and remained negligible until the subsequent overtopping. The cause of the difference between the field data and simulation data was not readily apparent. Again changes in the Ks values did not enable the model to replicate the behaviour of the field data. Both the field data and model simulation showed an increase in negative vertical head gradient (upward flow) as the shoreline approached with a peak at the point when the rising tide crossed

477

the monitoring location (Fig. 3). The field data showed a sharp change to positive gradient conditions (downward flow) immediately afterwards whereas the model predicted a more gradual change to a positive gradient. It is also apparent that the model was not able to accurately predict the depth to which the water table fell at both stations F and G, even with increased Ks values incorporated into the models. 3.1.3. Lower intertidal and subtidal bank edge – stations H, I and J The simulation data showed good agreement with the vertical head gradient data from monitoring stations (H, I and J) located in the lower intertidal and subtidal areas of the sandbank edge. At monitoring station H the vertical gradients predicted by the numerical model matched the dual-peak signal of negative

Fig. 4. Data for monitoring stations H and I (middle intertidal bank edge) comparing the observed head ((a) and (e)) and vertical head gradient ((c) and (g)) with the simulated head ((b) and (f)) and simulated vertical head gradient ((d) and (h)). All other presentation conventions are the same as in Fig. 2.

478

B. Gibbes et al. / Estuarine, Coastal and Shelf Science 80 (2008) 472–482

gradients (upward flow) measured at the field site (Fig. 4(a–d)). Both simulated and measured peaks in the negative gradient occurred when the shoreline crossed the monitoring point on the ebb and flood tide. The model was found to over-predict the magnitude of the negative head gradients (peak gradient z0.025 in field data cf. 0.045 in model data). Both the model and field data showed that flow was predominantly upward (negative head gradient) at this location indicating that it was a zone of net discharge of pore water across the sediment–water interface. Overall the simulated gradient at station I was consistent with the field data, both showing negative gradients (upward flow) throughout the exposure period (Fig. 4(e–h)). While the field data exhibited a double peak similar to those at station H, the model result only showed a single peak as the shoreline did not cross the monitoring point in the simulation. This was due to small differences between the low tide elevation at the field site and that of the synthetic tide used in the model. Station I was subject to a very short period of exposure in the simulation and hence resembled station J in the subtidal region. At station J both the field data and simulated data showed negative head gradients (upward flow) throughout the period when the upper parts of the sandbank were exposed (Fig. 5). This station was also characterized by a single peak of vertical head gradient that occurred at the low tide. The results indicated that this portion of the bank edge was also a discharge zone. At all three stations located on the lower parts of the sandbank edge (stations H, I and J in Fig. 4 of Gibbes et al. (2008)) the model appeared to slightly overestimate the peak negative gradients. The reason for this overestimation was not readily apparent but might be attributed to the inability of the model to replicate the development of a seepage face on the bank edge. 3.1.4. Summary of intra-tidal flows The above results focussed on the intra-tidal dynamics of the vertical head gradients. A key characteristic of such dynamics, as shown by both the field and simulation data, was that the vertical flow field was strongly centred on the shoreline. This was illustrated further in Fig. 6 where velocity vectors for the simulated pore water flow field were plotted at various points in the tidal cycle. The results showed that during the initial stages of sandbank exposure

on the ebb tide (Fig. 6(a)) the flow was directed predominantly towards the shoreline on the sandbank edge, although discharge towards the interior seagrass meadow was also evident. Around the low tide (Fig. 6(b)) the flow field was intensified below the bank edge and showed drainage of the whole sandbank with outflow occurring only around the shoreline on the bank edge. On the rising tide (Fig. 6(c)) the flow was again the strongest around the shoreline on the bank edge although a second flow system also existed near the intersection of the shoreline with the bank platform in the interior seagrass meadow. Near the point of complete overtopping on the rising tide (Fig. 6(d)) the flow remained the strongest near the shoreline location and a complex flow pattern was evident as the unsaturated area was filled. At all tidal stages the areas of strongest pore water flows were at the shoreline location. Simulated flow at the shoreline was predominantly vertical, consistent with the vertical head gradient data from the field measurements. The results shown in Fig. 6 also highlighted the temporal asymmetry of the flow field during the exposure period. For example, flow near the sandbank crest was predominantly downwards while upward flow dominated the lower intertidal zone of the bank edge. Such flow asymmetry will lead to residual pore water flow when averaged over the tidal cycle as discussed below. 3.2. Tide-averaged flow field To examine the tide-averaged flow field the mean hydraulic heads at points on a regularly spaced grid within the model domain were calculated by averaging the simulated head at each grid point over a single 12-h tidal cycle. The gradient of these averaged heads then determined the tide-averaged flow field. The result illustrates zones of infiltration and discharge across the sandbank (Fig. 7). The tide-averaged flow field near the bank edge was similar to that inferred from the measured head data of Gibbes et al. (2008) using the same averaging method. In particular, both the field and simulation data identified a net infiltration zone on the bank platform and a net discharge zone along the bank edge. Linking these two zones is a circulation system centred on the edge of the sandbank around the mean shoreline. Such a circulation system is a direct result (residual) of asymmetric intra-tidal flows discussed previously in Section 3.1.

Fig. 5. Data for monitoring station J (sub-tidal bank edge) comparing the observed head (a) and vertical head gradient (c) with the simulated head (b) and simulated vertical head gradient (d). All other presentation conventions are the same as in Fig. 2.

B. Gibbes et al. / Estuarine, Coastal and Shelf Science 80 (2008) 472–482

479

Fig. 6. Pore water flow field at various points during the tidal cycle: (a) early stages of bank exposure (elapsed time ¼ 129.1 h), (b) low tide (elapsed time ¼ 131.5 h), (c) rising tide (elapsed time ¼ 134.15 h) and (d) during overtopping (elapsed time ¼ 134.37 h). Arrows represent the instantaneous velocity vectors of the pore water flow field. The dashed line represents the elevation of the tidal water level. The water table location within the sandbank is represented by a solid thin line. The greyscale shading represents the effective saturation of the sandbank where a value of 1 represents completely saturated conditions and a value of 0 represents completely dry conditions.

The model predicted a second tide-averaged flow system below the bank platform (Fig. 7). In this interior flow system seawater is cycled from the bank platform to the interior seagrass meadow. Analysis of the simulated head gradient data for the interior seagrass meadow showed a similar pattern to that observed in the subtidal zone (station J, Fig. 5) on the sandbank edge (i.e., negative head gradients (upward flow) throughout the period when the upper parts of the sandbank were exposed). Although this interior circulation system was not measured in the field (monitoring equipment not placed in this part of the sandbank), its generation mechanism was the same as that for the first circulation system. Both circulation systems result from tide-induced asymmetric flows. The simulation results also showed a hydraulic divide at a horizontal distance of approximately 14 m, which separated these two flow systems. Results from the simulations that used the higher hydraulic conductivity values suggest that the location of this divide was influenced by the hydraulic conductivity. The divide location was found to shift away from the bank edge (i.e., moving towards the interior seagrass meadow) with increasing Ks values (divide location at horizontal distance z18 m and 23 m for Ks ¼ 4.05  104 m s1 and 8.101 104 m s1 respectively). It is speculated that this expansion of the bank edge circulation system is caused by increased drainage with increasing Ks values. With higher Ks values more water is able to discharge through the bank edge during sandbank exposure. This results in an expansion of the zone that feeds this discharge and a corresponding movement in the flow divide location occurs. This behaviour suggests that the bank edge flow system is dominant over the interior flow system. Analysis of the net discharge across the sediment–water interface during a tidal cycle (Fig. 8) also suggested that the bank edge flow system was more significant than the interior flow system in terms of pore water exchange with the Bay. The net discharge across the sediment–water interface in both flow systems were calculated by summing the fluxes across the sediment–water interface at 2 min intervals over a single 12-h tidal cycle. The result (Fig. 8) showed the influence of the topography of the upper

sandbank surface on the net infiltration with substantial variations across the bank platform. The predicted volume of pore water discharge across the sediment–water interface (per unit width) from the bank edge flow system over a tidal cycle (z 0.813 m3 m1 cycle1) was more than five times greater than that across the sediment–water interface in the interior seagrass meadow (z 0.143 m3 m1 cycle1). Transit times in the interior circulation system ranged between 2000 and 24,000 d with an average of approximately 11,000 d. These times are significantly longer than those in the bank edge circulation system which ranged from 30 to 2500 d with an average of approximately 182 d. The Darcian flow velocities within the bank edge circulation system were found to vary between 0.03 and 0.177 m d1 with a mean of approximately 0.092 m d1. This is similar to the range reported in many submarine groundwater discharge investigations

Fig. 7. Tide-averaged flow field predicted by the numerical simulations. The solid black lines represent the domain boundary while the dashed lines show the elevations of the high and low tide. The greyscale map represents the residence time of water entering through the upper bank platform. The solid white lines represent streamlines along the average flow path. Arrows represent the direction of the mean (residual) pore water flow.

480

B. Gibbes et al. / Estuarine, Coastal and Shelf Science 80 (2008) 472–482

Fig. 8. Simulated (a) tide-averaged specific pore water flux across the sediment–water interface (positive values indicate flow out of the sandbank). (b) low tide water table elevation/exit point location and (c) the slope of the sediment–water interface. In (b) the solid black lines represent the domain boundary while the dashed lines show the water levels at the high and low tide.

(0.05–1.0 m d1) (Burnett et al., 2006). In contrast the interior circulation system exhibited smaller Darcian flow velocities in the range of 0.005–0.015 m d1 with a mean of approximately 0.011 m d1. These rates are smaller than 0.05 m d1, the rate which is considered low or even marginally detectable in measurements of submarine groundwater discharge (Burnett et al., 2006). The transit times and Darcian pore water flow velocities were found to vary between models that used different Ks values. 4. Discussion 4.1. Simulation of field data The results show that while the numerical model was able to replicate the broad trends observed in the field data there were a number of specific features that were not able to be reproduced. In particular the model was unable to completely simulate surface water ponding on the upper sandbank platform as well as the extent of water table lowering. These differences could have implications in terms of accurately quantifying pore water exchange rates. Potential causes of these differences are discussed here with a view to identifying improvements to the modelling approach. The inability of the model to simulate the effects of microtopography on the drainage and trapping (ponding) of free surface water on the surface of the upper bank platform was evident in the vertical gradient data for the site in the middle of the bank platform (station C, Fig. 2). To simulate these micro-scale effects would require a better approximation of the drainage of free surface water on the upper sandbank surface. Accurate representation of surface water flows (particularly where the free-water surface is considered) and the interaction of these flows with porous medium is a challenging area for numerical modelling. Future model development to account for these processes would be a good direction

for additional research. Such an approach could be achieved via expanding the technique for sequentially coupling surface water and porous medium flow presented by Cardenas and Wilson (2006) to include losses from the surface water domain due to infiltration. The inability of the model to simulate the depth of water table lowering observed at the field site during sandbank drainage was assessed in a preliminary way via sensitivity testing based on three different Ks values. The simulation results suggests that, while the model is sensitive to Ks values, there are likely to be a range of other parameters such as soil moisture parameters (e.g., S, kr, C and Se) to which the model is also sensitive. However, given the relative homogeneity of the porous medium (Gibbes et al., 2008) and the relatively small range of variation in soil moisture properties for a sand sediment type (Carsel and Parrish, 1988), it is speculated that small-scale (0.05–0.10 m) topographic features are likely to exert a greater influence on pore water flow patterns at the site. This is supported by the simulation results which clearly show the area of greatest pore water flow is centred around the shoreline location, which in turn varies with tidal height and local topography. This suggests that sandbank geometry might exert a controlling influence on flow patterns. The more profound water table lowering observed in the field might also be due to three-dimensional drainage of the sandbank. In this case the two-dimensional numerical model would restrict the drainage of the sandbank resulting in a less pronounced lowering of the water table. Another source of variation might result from the inability of the model to accurately simulate seepage face dynamics (as discussed in Section 2.1.1). As the model only allowed flow across submerged portions of the sediment–water interface, drainage from the upper parts of the bank was forced to discharge further down the bank edge compared to a situation where a seepage face existed above the shoreline. Such forced discharge might be responsible for the model’s under-prediction of water table lowering in the bank platform as well as the over-prediction of the vertical head

B. Gibbes et al. / Estuarine, Coastal and Shelf Science 80 (2008) 472–482

gradients and hence flow rates in the lower intertidal and subtidal parts of the bank edge. Attempts to modify the boundary condition of the model to allow seepage face formation were unsuccessful (due to the introduction of numerical instabilities and added computation expense). Further model development is required to improve the predictive capacity of the model simulations. In particular further investigation of the relative influence of the porous medium properties (including the introduction of heterogeneities) compared to small topographic changes and three-dimensional flows is needed to better understand the importance of these features for future model development. As a first step an expanded sensitivity analysis covering a range of porous medium properties (e.g., Ks, S, kr, C, Se) as well as small random perturbations (0.05–0.10 m scale) in sandbank topography would assist in this regard. Accurate simulation of seepage face dynamics is also likely to improve the flow simulation. Notwithstanding the above limitations the current model has been useful in identifying broad scale pore water exchange patterns. 4.2. Pore water exchange patterns The simulation results suggested pore water exchange driven by the bank edge circulation system is concentrated in an area approximately 40 m wide on the sandbank edge (horizontal distance 0–40 m in Fig. 8). Such exchange is likely to result in export of pore water with elevated concentrations (compared to seawater concentration) of dissolved iron to the sediment–water interface in the discharge zone. In comparison, pore water exchange in the interior circulation system is relatively weak. Field data presented by Gibbes et al. (2008) showed a sharp peak of dissolved iron 0.4–0.6 m below the sediment–water interface at monitoring stations within the interior seagrass meadow (refer to data for stations A and B in Gibbes et al. (2008)). No measurable dissolved iron was detected near the sediment–water interface at these monitoring stations. The tide-averaged flow field predicted by the numerical model suggests that the interior seagrass meadow should be a zone of net pore water discharge. Under the influence of the interior circulation system, measurable dissolved iron concentrations would be expected near the sediment–water interface in this area, particularly when significant concentrations exist deeper in the sediment. However, the absence of measurable dissolved iron concentrations near the sediment–water interface suggests that other more important processes than the tidally induced flow are operating at this location. Near the sediment–water interface in the interior seagrass meadow, other pore water exchange mechanisms, such as currentand wave-induced subtidal pumping, could be more important. Darcian flow velocities due to current-induced advective exchange range from 0.24 to 800 m d1 (Thibodeaux and Boyle, 1987; Huettel and Gust, 1992; Huettel et al., 1996, 1998) whereas wave-induced advective exchange can range from 2.16 to 30 m d1 (Riedl et al., 1972; Webb and Theodor, 1972; Precht and Huettel, 2004). The upper bounds of these velocities are significantly higher than those predicted by the present tidal exchange model. Although subtidal pumping processes only operate when the sandbank surface is submerged, they are likely to complicate the interior circulation system identified here. Other exchange processes such as density-driven convective overturn due to temperature gradients (Rocha, 2000) could also be important. Pore water exchange due to density-driven convection occurs when permeable sediments that have been heated by solar radiation (during exposure) are flooded by relatively cold water on the rising tide. If this process occurs at the site, it is likely to operate for relatively short periods during the overtopping event. Further investigations are required to establish the significance of these potential pore water exchange processes at the field site.

481

The relatively long residence times of pore waters in both circulation systems, combined with the anoxic (reducing) condition in the sediments, allow for interactions between pore waters and ubiquitous solid phase forms of iron in the sediments. These interactions are likely to increase dissolved iron concentrations above those of the infiltrating oxygenated seawater. It is speculated that such an increase in dissolved iron concentrations might constitute an important local iron source for L. majuscula growth. Additional investigation into the coupled flow and reaction processes are required before specific conclusions can be made in this regard. 5. Conclusions and implications A two-dimensional numerical model was developed to simulate tide-induced pore water flow in an offshore intertidal sandbank in Moreton Bay, Australia. The simulation results generally agreed with the field data presented by Gibbes et al. (2008) for the same site. The simulated pore water flow systems presented here support the results of the field investigation by Gibbes et al. (2008) which suggest that tidally driven exchange might constitute a significant mechanism for transport of dissolved chemicals across the sediment–water interface. Therefore the current conceptual models of iron delivery to L. majuscula bloom sites in Moreton Bay (Dennison and Abal, 1999; O’Neil and Dennison, 2005) should incorporate tidally driven exchange as a potential delivery mechanism. The model also predicted an interior circulation system with net infiltration on the bank platform and net discharge to the interior seagrass meadow. This interior circulation system was not revealed by Gibbes et al.’s (2008) field investigation. Pore water exchange in this system is relatively weak compared with the flow system near the bank edge. In reality other pore water exchange processes, such as subtidal pumping, might have a greater influence on pore water flows in the interior seagrass meadow, particularly near the sediment–water interface. Differences between the simulation results and field data were observed. In particular the numerical models were found to underpredict the degree of water table lowering during sandbank drainage on the ebb tide. The model was also found to slightly overpredict the peak negative head gradients (upward flow) in the lower intertidal and subtidal areas of the bank edge. The reasons for these differences are unclear however it is suggested that the use of three-dimensional flow models, that are informed by detailed topographic information (cm scale) and that are capable of simulating seepage face dynamics, would offer an improvement over the approach used in this study. To the author’s knowledge this study is a first of its kind in combining modelling and field investigations on pore water flows in offshore intertidal sandbank systems. The findings presented here have implications for studies of not only offshore intertidal sandbanks but also nearshore salt marsh environments. The differences between the field data and numerical simulations highlight the complexity of the tidally driven exchange processes and the need for ongoing development and refinement of approaches to data collection and modelling for these intertidal systems. Acknowledgements This research was supported by an Australian Research Council Discovery Grant (DP0772660) and a National Natural Science Foundation of China Grant (no.: 50425926). The comments and suggestions of three anonymous reviewers are acknowledged. The financial assistance of a University of Queensland Graduate School Scholarship to the first author and an Australian Postgraduate Award to the second author is also gratefully acknowledged.

482

B. Gibbes et al. / Estuarine, Coastal and Shelf Science 80 (2008) 472–482

References Ahern, K.S., Ahern, C., Udy, J., 2007. Nutrient additions generate prolific growth of Lyngbya majuscula (cyanobacteria) in field and bioassay experiments. Harmful Algae 6, 134–151. Ahern, K.S., O’Neil, J.M., Udy, J.W., Albert, S., 2006. Effects of iron additions on filament growth and productivity of the cyanobacterium Lyngbya majuscula. Marine and Freshwater Research 57, 167–176. Albert, S., O’Neil, J.M., Udy, J.W., Ahern, K.S., O’Sullivan, C.M., Dennison, W.C., 2005. Blooms of the cyanobacterium Lyngbya majuscula in coastal Queensland, Australia: disparate sites, common factors. Marine Pollution Bulletin 51, 428–437. Anderson, P.A., Hunt, R.J., Krohelski, J.T., Chung, K., 2002. Using high hydraulic conductivity nodes to simulate seepage lakes. Ground Water 40, 117–122. Billerbeck, M., Werner, U., Bosselmann, K., Walpersdorf, E., Huettel, M., 2006. Nutrient release from an exposed intertidal sand flat. Marine Ecology Progress Series 316, 35–51. Burnett, W.C., Aggarwal, A., Aureli, A., Bokuniewicz, H., Cable, J.E., Charette, M.A., Knotar, E., Krupa, S., Kulkarni, K.M., Loveless, A., Moore, W.S., Oberdorfer, J., Oliveria, J., Ozyurt, N., Povinec, P., Privitera, A.M.G., Rajar, R., Ramessur, R.T., Scholten, J., Stieglitz, T., Taniguchi, M., Turner, J.V., 2006. Quantifying submarine groundwater discharge in the coastal zone via multiple methods. Science of the Total Environment 37, 498–543. Cardenas, M.B., Wilson, J.L., 2006. The influence of ambient groundwater discharge on exchange zones induced by current-bedform interactions. Journal of Hydrology 331, 103–109. Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water retention characteristics. Water Resources Research 24, 755–769. Celia, M.A., Bouloutas, E.T., Zarba, R.L., 1990. A general mass-conservative numerical solution for the unsaturated flow equation. Water Resources Research 26 (7), 1483–1496. COMSOL, 2004. FEMLAB 3.1 Multiphysics Modeling – User’s Guide. COMSOL Inc, Burlington, MA, USA, 305 pp. Dennison, W.C., Abal, E.G., 1999. Moreton Bay Study: a scientific basis for the healthy waterways campaign. South East Queensland Regional Water Quality Management Strategy, Brisbane, 246 pp. Gardner, L.R., 2005. A modeling study of the dynamics of pore water seepage from intertidal marsh sediments. Estuarine, Coastal and Shelf Science 62, 691–698. Gardner, L.R., Wilson, A.M., 2006. Comparison of four numerical models for simulating seepage from salt marsh sediments. Estuarine, Coastal and Shelf Science 69, 427–437. Gibbes, B., Robinson, C., Li, L., Lockington, D.A., Carey, H., 2008. Tidally driven pore water exchange within offshore intertidal sandbanks: Part I Field measurements. Estuarine, Coastal and Shelf Science 79, 121–132. van Genuchten, M.T., 1980. A closed form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal 44, 892–898.

Harvey, J.W., Germann, P.F., Odum, W.E., 1987. Geomorphological control of subsurface hydrology in the Creekbank Zone of Tidal Marshes. Estuarine, Coastal and Shelf Science 25, 677–691. Huettel, M., Gust, G., 1992. Impact of bioroughness on interfacial solute exchange in permeable sediments. Marine Ecology Progress Series 89, 253–267. Huettel, M., Ziebis, S., Forster, S., 1996. Flow-induced uptake of particulate matter in permeable sediments. Limnology and Oceanography 41, 309–322. Huettel, M., Ziebis, S., Forster, S., Luther, G.W., 1998. Advective transport affecting metal and nutrient distributions and interfacial fluxes in permeable sediments. Geochimica et Cosmochimica Acta 62, 613–631. Li, H., Li, L., Lockington, D.A., 2005. Aeration for plant root respiration in a tidal marsh. Water Resources Research 41, 1–11. W06023. Novakowski, K.I., Torres, R., Gardner, L.R., Voulgaris, G., 2004. Geomorphic analysis of tidal creek networks. Water Resources Research 40 W05401. O’Neil, J.M., Dennison, W.C., 2005. Chapter 6: Lyngbya. In: Abal, E.G., Bunn, S.E., Dennison, W.C., Healthy Waterways Healthy Catchments: Making the connection in South East Queensland, Australia. Moreton Bay Waterways and Catchments Partnership, Brisbane, pp. 119–148. Precht, E., Huettel, M., 2004. Rapid wave-driven advective pore water exchange in a permeable coastal sediment. Journal of Sea Research 51, 93–107. Reeves, H.W., Thibodeau, P.M., Underwood, R.G., Gardner, L.R., 2000. Incorporation of total stress changes in the ground water model SUTRA. Ground Water 38, 89–98. Riedl, R., Huang, N., Machan, R., 1972. The subtidal pump: a mechanism of intertidal water exchange by wave action. Marine Biology 13, 210–221. Robinson, C., Li, L., Barry, D.A., 2007. Effect of tidal forcing on a subterranean estuary. Advances in Water Resources 30, 851–865. Rocha, C., 2000. Density-driven convection during flooding of warm, permeable intertidal sediments: the ecological importance of the convective turnover pump. Journal of Sea Research 43, 1–14. Roy, P.S., Williams, R.J., Jones, A.R., Yassini, I., Gibbs, P.J., Coates, B., West, R.J., Scanes, P.R., Hudson, J.P., Nichol, S., 2001. Structure and function of South-east Australian Estuaries. Estuarine, Coastal and Shelf Science 53, 351–384. Thibodeaux, J.L., Boyle, J.O., 1987. Bed-form generated convective transport in bottom sediment. Nature 325, 341–343. Ursino, N., Silvestri, S., Marani, M., 2004. Subsurface flow and vegetation patterns in tidal environments. Water Resources Research 40, 1–11. W05115. Watkinson, A.J., O’Neil, J.M., Dennison, W.C., 2005. Ecophysiology of the marine cyanobacterium, Lyngbya majuscula (Oscillatoriaceae) in Moreton Bay, Australia. Harmful Algae 4, 697–715. Webb, J.E., Theodor, J.L., 1972. Wave-induced circulation in submerged sands. Journal of the Marine Biological Association of the United Kingdom 52, 903–914. Wilson, A.M., Gardner, L.R., 2006. Tidally driven groundwater flow and solute exchange in a marsh: Numerical simulations. Water Resources Research 42, 1–9. W01405. Winter, T.C., 1976. Numerical simulation and analysis of the interaction of lakes and groundwater. Transactions of the American Geophysical Union 57, 45.