1 Introduction

It is known that natural upwelling of deep water in the sea surface forms good fishery areas. The reason for this is that the enhanced vertical mixing by the upwelling transports nutrients from the bottom to the surface efficiently. As a result, primary production is activated in the surface and, consequently, a good fishery can be expected. An artificial mound settled at the sea bottom possibly produces similar upwelling artificially. This technology is expected to contribute not only to the fishery, but also to the biological fixation of CO2 from the atmosphere to the ocean. Compared with other engineering methods to enhance vertical mixing, such as water pumping or air bubbling, an artificial mound has an advantage in a sense that it does not consume external energy because it converts energies of tides or currents into vertical mixing.

It was reported that an artificial structure settled at the sea bottom alters flow pattern and surrounding ecosystems (e.g., [13]). It was also suggested from field observation that natural seamounts enhance vertical mixing. Polzin et al. [7], Lerwell et al. [2], and St. Laurent et al. [11] reported large vertical diffusivity above rough topography in the deep ocean. Lueck and Mudge [3] found large mixing near a shallow seamount at most 104 times as much as that of the site far away from it. Toole et al. [12] found that the vertical diffusivity is enhanced within 500 m vertically and 1–3 km horizontally from the flank of the seamount. The effect of a seamount on vertical mixing depends on both the environmental conditions, such as currents and stratification, and the specification of a mound and numerical approach is useful to quantify it.

In this study, we focused on the physical process of the mixing effect by an artificial mound settled in the shallow ocean and reproduce the turbulent flow field around the mound by a numerical simulation method. It is considered that the mechanism of the enhancement in vertical mixing is attributed to breaking of internal waves generated by the interaction between tidal currents and the mound. In this process, turbulence plays a significant role and the large eddy simulation (LES) was applied to express its effect. Energy dissipation rate was calculated from the simulation results and compared with observation in the real sea. Vertical diffusivity was also estimated, and its relation with current velocity and stratification was investigated.

2 Measurement of velocity and stratification at artificial mound

2.1 Artificial mound and flow measurement

An artificial mound was settled at 33°28′41″N and 129°25′29″E, north-northwest to Ikutsuki Island (Fig. 1). The depth of water at this sea site is about 80 m, and the direction of the main tide is about 23° to 203°. The mound consists of mass of rectangular concrete blocks made of coal ash. Figure 2 shows schematic top view of the mound settled at the sea bottom. Sonar observation revealed that the mound elongates in the south-southeast and north-northwest direction with two peaks. The height of each top is about 12 m, and the horizontal length is 60 m for each crest.

Fig. 1
figure 1

Bathymetry around the artificial mound. The circle indicates the place where the mound was set, 33°28′41″N and 129°25′29″E

Full size image
Fig. 2
figure 2

Schematic diagram of the location of the mound and the fixed ADCP used in the measurement, a top view and b side view. y indicates the cross-flow direction

Full size image

From July 28 to September 13, 2005, vertical profiles of flow velocity were measured by using an ADCP (RDI Sentinel, 300 kHz) set at an edge of the mound. The spatial resolution of the vertical profile was 2 m, and the sampling period was 180 s. Stratification was also measured on 18 August and 13 September by a CTD.

Figure 3 shows time series of velocity in the main tidal direction. About 2-week spring–neap tidal cycle was observed. Spring and half tides correspond to 18 August and 13 September, respectively. Hereafter, we call the former and the latter as Case FS and Case FH, respectively, for distinction. The amplitudes of velocity were approximately 0.8 and 0.4 m s−1 in Cases FS and FH, respectively. Figure 4 shows the vertical profiles of density and buoyancy frequency obtained from the CTD measurement. The stratification intensity in Case FS is found to be larger than that in Case FH by a factor of 2.

Fig. 3
figure 3

Horizontal velocity in the tidal direction obtained by the fixed ADCP

Full size image
Fig. 4
figure 4

Ensemble-averaged vertical profiles of a density and b buoyancy frequency in Cases FH and FS

Full size image

2.2 Estimation of energy dissipation rate from velocity profiles

Energy dissipation rates were estimated from the time series of velocity profiles measured by use of the ADCP fixed to the sea bottom. The spatial interval of velocity profile was 2 m, which is not sufficiently small for the direct estimation of turbulence dissipation in the way done by a microstructure profiler. Alternatively, we adopted the indirect estimation model of Mackinnon and Gregg [4] (hereafter, the MG model). The MG model is a modification of the model of Gregg [1] (hereafter, the GH model) to apply to shallow coastal seas and represents the energy dissipation rate as

$$ \varepsilon_{\text{MG}} = \varepsilon_{0} \left( {\frac{N}{{N_{0} }}} \right)\left( {\frac{{S_{\text{lf}} }}{{S_{0} }}} \right) $$
(1)

where ε0 is the reference energy dissipation rate which was set to be 6.9 × 10−10 m2 s−3 after Mackinnon and Gregg [4], N is the buoyancy frequency, and Slf is the low-frequency component of 4-m shear filtered by 0.17 cph (2.97 × 10−4 s−1). N0 and S0 are the reference buoyancy frequency and vertical shear, respectively, both of which are set to be 3 cph (5.24 × 10−4 s−1).

3 Numerical simulation of turbulence dissipation around artificial mound

3.1 Governing equations

The LES of flow field around an artificial mound was carried out. The governing equations are the continuity equation, Navier–Stokes equation with the Boussinesq approximation, and the density transport equation.

$$ \frac{{\partial u_{j} }}{{\partial x_{j} }} = 0 $$
(2)
$$ \frac{{\partial u_{i} }}{\partial t} + \frac{{\partial \left( {u_{i} u_{j} } \right)}}{{\partial x_{j} }} = - \frac{1}{{\rho_{0} }}\frac{\partial p}{{\partial x_{i} }} + \frac{\partial }{{\partial x_{j} }}\left[ {\left( {\nu_{0} + \nu_{\text{SGS}} } \right)\frac{{\partial u_{i} }}{{\partial x_{j} }}} \right] - \frac{{\rho^{\prime } g}}{{\rho_{0} }}\delta_{i3} + f_{i} $$
(3)
$$ \frac{{\partial \rho^{\prime } }}{\partial t} + \frac{{\partial \left( {\rho^{\prime } u_{j} } \right)}}{{\partial x_{j} }} = \frac{\partial }{{\partial x_{j} }}\left[ {\left( {\frac{{\nu_{0} }}{\text{Pr}} + \frac{{\nu_{\text{SGS}} }}{{\text{Pr}_{\text{SGS}} }}} \right)\frac{{\partial \rho^{\prime } }}{{\partial x_{j} }}} \right] + u_{3} \frac{{N^{2} \rho_{0} }}{g} $$
(4)

where ui is the velocity vector, p is the pressure, ν0 is the kinematic viscosity, νSGS is the subgrid-scale eddy viscosity, g is the gravity acceleration, ρ0 is the reference density, ρ′ is the density from which the background density stratification is subtracted and normalized by ρ0, f is the driving force, Pr is the Prandtl number, and PrSGS is the turbulent Prandtl number. The last term in the right-hand side of Eq. 4 is corresponding to the forcing of background stratification. A tidal acceleration was imposed as a driving force, f, by a cosine function with amplitude U0 and frequency ω.

$$ f_{i} = \delta_{i1} U_{0} \omega \cos \omega t $$
(5)

3.2 Turbulence model

The characteristic length and velocity scales in the present study are represented by the height of the mound and the amplitude of the tidal flow, respectively. The Reynolds number based on these scales becomes about 9.6 × 106 and 4.8 × 106 for spring and half tides, respectively. In the present LES, eddies smaller than the grid scale are expressed by an eddy viscosity model of Smagorinsky [9]:

$$ \nu_{\text{SGS}} = \left( {C_{S} \Delta } \right)^{2} \left( {2S_{ij} S_{ij} } \right)^{1/2} $$
(6)

where CS is the Smagorinsky constant which was set to be 0.1, Δ is the grid size, and Sij is the strain-rate tensor.

3.3 Numerical procedure

Equations 2–4 were solved in a time-developing way by a finite volume method. The computational domain used was divided by the body-fitted grid system. The velocity and pressure were allocated at the cell center, and the method of Rhie and Chow [8] was adopted to avoid the checker-board instability. Spatial accuracy was fourth order for the convection term and second order for the other terms. A third-order Adams–Bashforth method was used for time accuracy.

Periodic boundary condition was imposed in the horizontal directions, and free-slip boundary condition was adopted at the top and bottom boundaries. It is considered that the kinetic energy of fluctuating components is transported by eddies and internal waves. Eddies are associated with the bottom boundary layer and the bottom topography, and the scale of the former is considerably smaller than that of the latter. Although the effect of eddies in the bottom boundary layer is large in the region near the sea bottom, its effective distance may be limited compared with those of topographical eddies and internal waves. It is deduced that the friction effect is little in the region far from the topography and was not considered in this study.

The length scales of computational domain were 60 m × 60 m in the horizontal directions and 80 m in the vertical direction. The mound is represented by a vertically axisymmetric cosine curve with a height and a wavenumber of 12 m and 1/60 m−1, respectively. Figure 5 shows the vertical cross section of computational domain. The number of total grids was 64 × 64 × 64 = 262,144. The spatial resolutions were 0.94 m in the horizontal directions and 0.24–3.60 m in the vertical direction.

Fig. 5
figure 5

Calculation domain and grids (side view)

Full size image

3.4 Calculation conditions

Table 1 shows the fundamental physical parameters used in the present simulation. Two cases of simulations are carried out under the two sets of conditions of U0 and N, the values of which are representative of those measured on August 18 and September 13, 2005, as shown in Table 2. In Case FS, N is given by a linear function which provides 1.8 × 10−2 and 3.0 × 10−2 s−1 at the sea bottom and the sea surface, respectively.

Table 1 Common parameters
Full size table
Table 2 Calculation cases
Full size table

Time period of reciprocal flow that mimics a single tidal constituent was assumed to be semidiurnal, T = 43,200 s.

Each case was calculated for the period corresponding to three tidal cycles, in which the first two cycles are acceleration stage. The last cycle was divided into four timings, namely t = T/4, T/2, 3T/4, and T. Theoretically, the velocity of tidal flow with no disturbance becomes its maximum at t = T/4 and 3T/4 and becomes 0 at t = T/2 and T. In the following analysis, results at t = T/4 and t = T/2 are mainly focused.

4 Results and discussions

4.1 Flow field

Figure 6 shows the vertical profiles of calculated horizontal velocity in the x direction in Cases FS and FH at t = T/4 and t = T/2. The velocity was averaged on the horizontal plane, and the uniform tidal flow was subtracted. It is seen that the velocity does not become 0 even at t = T/2 especially near the bottom. This means that the effect of bottom topography on the vertical shear is not negligible. Figure 7 shows the contour maps of calculated w on the vertical section across the center of the mound at t = T/4. Upward and downward flows can be found at the upstream and downstream sides of the mound, respectively. The amplitudes of upward and downward velocities in Case FS are almost twice as large as those in Case FH, because of the amplitude difference of forced tides. The large upward velocity can only be seen in the region <20 m, and it is suggested that the effect of topographical eddies is limited near the bottom topography.

Fig. 6
figure 6

Vertical profiles of streamwise velocity from which mean velocity, U, was subtracted, in Cases FS (a) and FH (b). Solid and dashed lines show profiles at t = T/4 and T/2, respectively

Full size image
Fig. 7
figure 7

Vertical distributions of vertical velocity at t = T/4, a Case FS and b Case FH

Full size image

Figure 8 shows the spatial variation of density fluctuation on the same cross section as Fig. 7 at t = T/4 in Cases FS and FH. Wavelike structures can be seen in the upstream of the mound in the both cases. Considering that the scale of the tidal excursion in a cycle is much larger than that of the mound, they are corresponding to the quasi-steady lee waves as classified by St. Laurent and Garrett [10]. The slope of the internal wave is estimated by ω2/(N2 – ω2), and it is definitely smaller than that of the mound. Therefore, the topography in this study is supercritical. It can be seen in Fig. 8 that the effective range of internal waves is 40–50 m and it is larger than that of topographical eddies.

Fig. 8
figure 8

Vertical distributions of density fluctuation at t = T/4, a Case FS and b Case FH

Full size image

4.2 Energy dissipation rate

Energy dissipation rate was calculated from the LES results by using

$$ \varepsilon = 2\left( {\nu_{0} + \nu_{\text{SGS}} } \right)S_{ij} S_{ij} $$
(7)

Figures 9 and 10 show the vertical distributions of energy dissipation rate in Cases FS and FH, respectively, at t = T/4 and T/2, at which the tidal velocity is at the maximum and 0, respectively. Large energy dissipation can be found above the mound at t = T/4 in the both cases. On the other hand, the similar trend cannot be seen at t = T/2. This suggests that the energy dissipation above the mound should be attributed to the shear produced by the internal waves shown in Fig. 8. The amplitude of the dissipation is correlated with the tidal speed, as can be seen by comparing Cases FS and FH.

Fig. 9
figure 9

Vertical distributions of energy dissipation rate in Case FS, at = T/4 and bt = T/2

Full size image
Fig. 10
figure 10

Vertical distributions of energy dissipation rate in Case FH, at = T/4 and bt = T/2

Full size image

Figures 11 and 12 show the horizontal distributions of energy dissipation rate at z = 20 m in Cases FS and FH, respectively, at t = T/4 and T/2. At t = T/4, large energy dissipation can be seen on the upstream side of the mound top and the spatial variation is large in the both x and y directions. At t = T/2, on the other hand, the spatial distribution is relatively small. These results suggest that the temporal variation of the energy dissipation rate is large on the upstream side of the mound top where large energy dissipation can be found at t = T/4, while it is small in the other region.

Fig. 11
figure 11

Horizontal distributions of energy dissipation rate at z = 20 m in Case FS at at = T/4 and bt = T/2

Full size image
Fig. 12
figure 12

Horizontal distributions of energy dissipation rate at z = 20 m in Case FH at at = T/4 and bt = T/2

Full size image

4.3 Comparison of energy dissipation rate with measurement

Figure 13 shows the vertical profiles of energy dissipation rate obtained from the LES results. These were obtained by taking the spatial average in the horizontal plane over the computational domain. Both the GH and MG models using measured data are also plotted in Fig. 13. It is found that the GH model is always larger than the MG model, and this agrees with the discussion of Mackinnon and Gregg [4] that the GH model overestimates the energy dissipation rate in the shallow ocean. It is noted that the energy dissipation rate is larger in Case FS than Case FH for the empirical MG model as well as LES results. The effect of tidal speed on the energy dissipation rate is significant. The orders of magnitude of the LES result agree with that by the MG model at t = T/2, at which tidal speed is almost 0, but is larger at t = T/4, at which tidal speed is its maximum. Considering that the estimation by the MG model in Fig. 13 is the temporal mean, the present LES tends to overestimate energy dissipation rate when tidal speed is nonzero.

Fig. 13
figure 13

Vertical distributions of energy dissipation rate obtained from LES results and observation, a Case FS and b Case FH. Solid and dashed lines show profiles at t = T/4 and t = T/2, respectively. Open circles and triangles show estimates by GH and MG models, respectively

Full size image

A possible reason for this discrepancy is that the location of the ADCP may not be within the region of large energy dissipation. As can be seen in Figs. 11 and 12, the energy dissipation rate has a large spatial variation at t = T/4 and its spatial mean on the horizontal plane is influenced by the region of large energy dissipation. Namely, the spatial mean tends to be an overestimation for the other region. At t = T/2, on the other hand, there is relatively small spatial variation and the spatial mean can represent the local energy dissipation in any region. As shown in Fig. 2, the location of the ADCP was off the tidal pass though the mound top, and the energy dissipation rate may be small with little temporal variation.

4.4 Vertical diffusivity

Osborn [6] suggested the vertical diffusivity of the water patch, the length scale of which is O (10 m), as

$$ K_{\rho } = \varGamma \frac{\varepsilon }{{N^{2} }} $$
(8)

where Γ is the mixing efficiency, usually set to be 0.2 as suggested by Osborn [6] and Moum [5]. Figure 14 shows the vertical profiles of Kρ calculated from the LES results by using Eq. 8. Kρ ranged from the order of 10−6–10−3 m2 s−1, depending on the vertical location and the tidal speed. There can be seen a large difference in Kρ between t = T/4 and t = T/2. Kρ at the former is more than ten times as large as that at the latter. It is considered that this difference is mainly resulted from the difference of energy dissipation rate.

Fig. 14
figure 14

Vertical profiles of vertical diffusivity in Cases FS and FH, calculated from the LES results

Full size image

It is noted that the diffusivity in Fig. 14 is the spatial mean on the horizontal plane, and its vertical distribution at t = T/4 is strongly influenced by the locally large energy dissipation near the mound top, as can be seen in Figs. 11 and 12. On the other hand, diffusivity at t = T/2 has little horizontal variation, ant it is less influenced by the mound. In other words, diffusivity at t = T/2 can be regarded as the background diffusivity without the effect of the mound.

4.5 Estimation of nutrients transport

Transport rate of nutrients from the bottom to the surface of the ocean can be quantified by using the vertical diffusive velocity. The diffusive velocity, wdiff, is given by

$$ w_{\text{diff}} = - K_{\rho } \frac{\partial C}{\partial z} $$
(9)

where C is the dimensionless concentration of nutrients. Although the actual transport depends on the local distributions of nutrients and diffusivity, it is roughly estimated using the vertical diffusivity estimated in the present simulation. For simplicity, it is assumed that the nutrient is distributed linearly in the vertical direction with dimensionless concentration being 1 at the sea bottom and 0 at the sea surface. At t = T/4 when the tidal speed is its maximum, Kρ becomes approximately 10−3 m2 s−1 near the sea bottom and 10−4 m2 s−1 near the sea surface. Therefore, wdiff results in the order of 10−5–10−6 m s−1 near the bottom and the surface, respectively. At t = T/2, when the tidal speed is almost 0, on the other hand, Kρ is approximately 10−5 m2 s−1 and wdiff results in 10−7 m s−1. Assuming that Kρ at t = T/2 is corresponding to the background diffusivity with little effect of the mound, a simple dimensional analysis shows that the time required for nutrients transportation from the bottom to the surface is 103 days. It can be said that installation of the mound accelerates this process by more than ten times faster.

5 Conclusion

In this study, flow fields around an artificial mound, the height of which is 12 m, in the shallow ocean, the depth of which is 80 m, were numerically simulated by an LES method, and its results were compared with the measurement. We focused on two cases with different tidal amplitude and stratification. Large energy dissipation attributed to topographical eddies can be seen near the mound, but its effect is limited in the range of about 20 m above the mound. It is implied that the energy dissipation in the range far from the mound is caused by the breaking of internal waves. The energy dissipation was large in the case of spring tide, in which both tidal amplitude and stratification are large. The energy dissipation rate estimated from LES results has large spatial variation when the tidal velocity is its maximum, while little spatial variation can be seen when the tidal velocity is approximately 0. The vertical profile of the latter agrees with that by the measurement, and it is deduced that the location of the ADCP was out of the region of large energy dissipation enhanced by the mound. The vertical diffusivity estimated from the LES results has a variation by a factor of 10 depending on the tidal speed, which may be regarded as the mound effect itself.

The present simulation can be used to predict the amount of nutrients transported to the ocean surface and how much primary production and CO2 fixation are expected. The rough estimation shows that the background diffusive velocity of nutrients was 10−7 m s−1 and it becomes more than ten times larger with the mound effect. Further prediction for primary production and consequent CO2 fixation requires consideration of chemical and biological processes which will be another work. The present simulation is also useful for designing the optimum mound shape under certain sea sites. The numerical model used in the present simulation can be applied to both the shallow and the deep ocean. It should be noted that, as the water depth becomes shallower, the top and the bottom boundary effects will become significant and application of the more relevant boundary conditions such as moving boundary and frictional boundary may be necessary.