Upwelling of deep seawater to the region, where sunlight reaches, can produce the ocean farm since deep seawater contains high concentration of nutrient. The numerical simulation for upwelling of deep seawater with the perpetual salt fountain proposed by Stommel et al. was conducted in this study. The temperature and salinity distributions measured in Mariana area where the upwelling experiment was conducted by Maruyama et al. was used. As a result, the velocity profile of the upwelling experiment was predicted as M-shape flow and the flow rate was estimated as 43t/day in the pipe. Additionally the possibility of reverse flow in the pipe was indicated. Furthermore the possibility of upwelling in other ocean areas using the results was discussed. As a result, it became clear that the unified representation of ocean conditions was achieved by the new dimensionless number RaR, which was modified Rayleigh number, and flow rate in the pipe could be evaluated by RaR

Komiya, A.

Maruyama, S.

Sato, T.

Tsubaki, K.

University of Queensland eSpace

16th Australasian Fluid Mechanics Conference 
Crown Plaza, Gold Coast, Australia 
2-7 December 2007 
 
Numerical Simulation of Upwelling Flow in Pipe Generated by 
Perpetual Salt Fountain 
 
Tetsuya SATO1, Shigenao MARUYAMA2, Atsuki KOMIYA2 and Koutaro TSUBAKI3 
1Graduate School of Engineering 
Tohoku University, Sendai, Miyagi, 980-8579 JAPAN 
2Institute of Fluid Science 
Tohoku University, Sendai, Miyagi, 980-8577 JAPAN 
3Department of Mechanical Engineering 
Saga University, Saga, Saga, 840-8502 JAPAN 
 
Abstract 
Upwelling of deep seawater to the region, where sunlight reaches, 
can produce the ocean farm since deep seawater contains high 
concentration of nutrient. The numerical simulation for upwelling 
of deep seawater with the perpetual salt fountain proposed by 
Stommel et al. was conducted in this study. The temperature and 
salinity distributions measured in Mariana area where the 
upwelling experiment was conducted by Maruyama et al. was 
used. As a result, the velocity profile of the upwelling experiment 
was predicted as M-shape flow and the flow rate was estimated 
as 43t/day in the pipe. Additionally the possibility of reverse flow 
in the pipe was indicated. Furthermore the possibility of 
upwelling in other ocean areas using the results was discussed. 
As a result, it became clear that the unified representation of 
ocean conditions was achieved by the new dimensionless number 
RaR, which was modified Rayleigh number, and flow rate in the 
pipe could be evaluated by RaR. 
 
Introduction  
The great increase of the world population has required greater 
food productivity. However, an expansion of farm production has 
been limited by the deforestation problem. In contrast, pelagic 
zone has large areas but has low biological productivity. These 
are classified as ocean desert. The fertilization of pelagic zone 
may solve a problem with food product. 
 
One of the methods for the fertilization of ocean desert is an 
effective utilization of deep seawater. Deep seawater, which is 
deeper than 200m, has rich nutrient concentration as compared to 
epipelagic seawater [6]. Therefore if it is possible to draw up 
deep seawater to photic zone, the sea cultivation can be created in 
ocean desert. 
 
Generally a pump has been used for upwelling of deep seawater. 
For example, the upwelling experiment has been conducted in 
Sagami bay, Japan [7]. However a pump requires a large amount 
of energy and a periodic maintenance. Furthermore the drawn up 
water may sink because it is drawn up with keeping low 
temperature. Therefore a pump is not suitable system for the 
artificial upwelling in pelagic zone. 
 
Maruyama et al. proposed the upwelling of deep seawater in 
pelagic zone using the perpetual salt fountain proposed by 
Stommel et al. [1,5]. In many tropical and subtropical regions, 
salinity and temperature of epipelagic seawater are higher than 
those of deep seawater. When a pipe is settled in such 
stratification and then filled with low salinity deep seawater, the 
inside temperature distribution along the pipe axis becomes same 
as outside one eventually. As a result, buoyancy force occurs in 
the pipe because inside salinity is lower than outside one. The  
 
 
Figure 1. Concept of the perpetual salt fountain 
 
concept of the perpetual salt fountain is shown in Figure 1. This 
buoyancy force induces the upwelling flow in the pipe and it 
continues as long as temperature and salinity difference between 
epipelagic and deep seawater exist. The artificial upwelling does 
not need external energy except the initial input. 
 
Tsubaki et al. conducted upwelling experiment in 2002 using the 
perpetual salt fountain in Mariana area, which lies in 142.24E 
and 11.25N of pelagic zone. A tracer was ejected in the pipe and 
one tracer sensor was laid 2m downstream of ejector to observe 
the upwelling flow. The upwelling velocity and tracer diffusivity 
was measured for the first time in the world [5]. However 
measured diffusivity was 104 times as large as molecular 
diffusivity. They predicted that this was the effect of many 
vortices produced by deformed and oscillated pipe.  
 
Zhang et al. conducted numerical simulation related to ocean 
experiment to predict flow field in the pipe [8]. The measured 
diffusivity was specified as eddy diffusivity in that calculation. 
As a result, agreement between numerical and experimental 
results was observed and velocity profile in the pipe was 
estimated as M-shape flow. However they compared the 
numerical result with the experimental one at only one point 
because one sensor was installed in the experiment. Therefore 
numerical result of velocity profile was partially reliable. 
 
Accordingly in this paper, the numerical simulation was 
conducted again and the results were compared with 
experimental one at two sensor points, which were laid 2m and 
5m downstream of the ejection. Then exact flow rate estimation 
at Mariana area was discussed. These sensor data were obtained 
after foregoing calculation [2]. Furthermore, as a first step for 
estimation of sea cultivation in other ocean area, the possibility of 
flow rate estimation in other ocean area was studied.  
394
Numerical Simulation 
Figure 2 is the calculation domain considered in this study. This 
domain shows a portion of ocean including the pipe and the top 
boundary represents ocean surface. The left boundary is specified 
as the pipe axis and the top and bottom boundary are specified as 
adiabatic wall. The right boundary is specified as the wall, which 
has ocean temperature. The pipe length and location correspond 
to ocean experiment. This domain is set large enough to take no 
notice of wall effect. 
 
 
 
Figure 2. Calculation domain and coordinate system 
 
The governing equations are as follows. 
( ) 0ρ∇ ⋅ =V ,                   (1) 
2( ) ( ) P
t
ρ ρ ν ρ∂ +∇ ⋅ = ∇ −∇ +∂
V VV V g ,     (2) 
2( ) ( )
p
T kT T
t c
ρ ρ∂ +∇⋅ = ∇∂ V
,           (3) 
    ( ) ( ) ( )C C D C
t
ρ ρ ρ∂ +∇⋅ = ∇⋅ ∇∂ V ,          (4) 
where , , , pk cρ ν  are density, kinematic viscosity, thermal 
conductivity and specific heat, respectively. The vectors V, g are 
velocity and gravity, T, C, D are temperature, concentration and 
mass diffusivity, respectively. 
 
The fluid outside the pipe is specified as steady temperature and 
salinity distributions, which were measured in ocean experiment 
in 2002, as shown in Figure 3 [5]. The fluid inside the pipe is 
specified as constant salinity, which was measured at the bottom 
of the pipe. At initial time, temperature distribution inside the 
pipe is specified as that outside the pipe. 
 
In this study, the measured diffusivity 10-5m2/s is assumed to be 
eddy diffusivity. Therefore eddy diffusivity becomes much larger 
than molecular diffusivity. As a result, eddy thermal diffusivity 
and eddy kinetic viscosity can be assumed to be same as eddy 
diffusivity. In this study, the mass diffusivity, D, in eq. (4) is 
specified as the eddy diffusivity. Furthermore the thermal 
diffusivity and kinematic viscosity are also specified as the eddy 
diffusivity. Thermal conductivity and viscosity, which are 
calculated with eddy diffusivity, are used in eqs. (2) and (3). The 
density is calculated by UNESCO’s equation, which assumes 
density as a function of temperature, salinity and pressure [4]. 
Other properties are assumed to be constant.  
 
 
Figure 3. Temperature and salinity distributions in Mariana area 
 
The numerical simulation code is built with FLUENT 6.3.13, a 
commercial CFD package. An unsteady simulation is conducted 
with first order implicit method. A structured non-uniform grid is 
used. The convection terms are modelled with second-order 
upwind scheme and the diffusion terms are modelled with 
second-order central scheme. The SIMPLE algorithm is used for 
pressure correction.  
 
There are two processes to estimate the flow rate. The first one is 
the estimation of velocity profile in the pipe using the calculation 
described above. The second one is the calculation of time-series 
tracer diffusion using the new calculation domain, which focuses 
on only the inside pipe around sensor and ejection. The velocity 
profile estimated by first step is applied to that calculation 
domain and assumed to be stable. Then, comparison of time-
series tracer concentration between numerical results and 
experimental one, which was obtained in 2004 [2], is conducted. 
Re-examination of steady velocity profiles is conducted if 
disagreement is observed. 
 
Numerical Results 
Flow Rate Estimation in Mariana Area 
The calculation results of axial velocity profiles at several typical 
horizontal levels are shown in Figure 4. The velocity profiles are 
observed as M-shape flow except for distribution at the inlet of 
the pipe. This is because heat transfer from outside fluid does not 
reach to the pipe axis. The flow rate in Figure 4 is 166.8t/day. 
Figure 4 also shows that the reverse flow may occur near the axis 
at middle of the pipe. In the experiments, upwelling flow was 
measured continuously at the sensor point, which lied 13m below 
the pipe outlet. However downward flow might be measured if 
the sensors were laid at different point. Therefore axial velocity 
at the axis is investigated to confirm reverse flow area and shown 
in Figure 5. The reverse flow between 20m and 190m below the 
pipe outlet is shown.  
 
Then, flow rate estimation is conducted using the method 
described above. The velocity profile of sensor point shown in 
Figure 4 is used. As a result, the comparison of time-series tracer 
concentration between numerical and experimental result is done 
and disagreement between the numerical simulation and the 
experimental one is obtained. The peak of time-series tracer 
concentration by numerical simulation is observed much earlier 
than that by experiment at two sensor points as shown in Figure 6. 
This indicated that the calculation results of velocity profiled 
shown in Figure 4 are over estimated.  
395
 
 
Figure 4. Axial velocity profiles at several levels in the pipe 
 
 
 
      Figure 5. Axial velocity distribution at the axis of the pipe 
 
Two factors are considered as the reason of over-estimation. The 
first factor is over-estimation of eddy diffusivity in the pipe. In 
this study, thermal conductivity and viscosity are estimated using 
the constant eddy diffusivity. However in reality, laminar sub 
layer exists and eddy diffusivity becomes less near the wall. 
Therefore this calculation may over estimates the heat transfer 
from outside the pipe. The second factor is inclination of the pipe. 
Three depth indicators were installed at top, middle and bottom 
of the pipe in experiment. In 2004 experiment, the measured 
value of depth at middle and bottom of the pipe was fluctuating 
even though the top of the pipe was almost stable [2]. This 
indicates that the pipe inclined. Therefore sensor location might 
be moved from the axis of the pipe and measurement point of the 
numerical simulation may not be same as that of experiment.  
 
Accordingly for the overcome of first factor, eddy diffusivity is 
reconsidered. However it is difficult to provide distribution of 
eddy diffusivity, which includes wall effect, because there are no 
enough experimental data near the wall. Therefore velocity 
calculated above is reconsidered and velocity magnitude is 
modified to coincide with experimental data keeping M-shape 
flow. Furthermore for the overcome of second factor, the 
measurement points of tracer concentration are repositioned 
radially. Considering the depth indicator data and sensor length, 
the measurement point was predicted to move about 5cm far 
from the axis in 2004 experiment. Therefore new measurement 
points are positioned at 5cm, 6cm and 7cm far from the axis. As a 
result, it becomes clear that the time-series tracer concentration 
obtained by quarter of velocity profile calculated above is good 
agreement with experimental data at 7cm far from the axis as 
shown in Figure 7. This indicates that velocity profile in the 
experiment may be close to M-shape flow. Furthermore, the flow 
rate of 2004 experiment in Mariana area is predicted as 43t/day 
 
 
Figure 6. Comparison between measurement data and calculated data 
 
 
 
Figure 7. Comparison between measurement data and corrected 
calculated data under the sensor inclined condition 
 
by using modified velocity profile. This value is 4 times less than 
the case in Figure 4. This flow rate depends not only on 
temperature and salinity but also on wave conditions. Therefore it 
may decrease at same ocean area if the effect of pipe oscillation 
and deformation are small. Accordingly measurement diffusivity 
10-5m2/s is assumed as maximum eddy diffusivity here.   
 
Flow Rate Estimation in Other Ocean Area 
As described above, there is an urgent need to create sea planting 
in pelagic zone because of problem with food product recently. 
Therefore not only Mariana area but also other pelagic areas 
should be paid attentions. However upwelling flow rate by 
perpetual salt fountain depends on ocean conditions and it is 
uncertain that enough flow rate can be obtained in all pelagic 
zone. Accordingly rough flow rate should be predicted in several 
pelagic zones before experiment. Therefore this section focuses 
on the possibility of flow rate estimation depending on ocean 
conditions numerically and discusses it with paying attention to 
101.29E and 19.49S ocean area near Australia particularly. 
 
The governing equations and calculation scheme are same as 
Mariana area case. The temperature and salinity distributions of 
the ocean area, which are obtained by GEOSECS (Geochemical 
Ocean Sections Study), are used as initial and boundary 
conditions. The pipe length is modified as 600m because the 
depth difference between maximum and minimum salinity is 
600m. The depth of outlet pipe is same as that in Figure2. The 
assumptions are also same as previous one. In this calculation, 
molecular and measured diffusivity 10-5m2/s are used as 
minimum and maximum flow rate case. Several pipe diameters, 
20cm, 30cm, 50cm, 100cm, are considered for the comparison. 
396
The calculation result of flow rate at each condition is shown in 
Figure 8. The calculated flow rate increases with the increment of 
pipe diameter and increment rate of measured diffusivity is much 
larger than that of molecular diffusivity. This is the effect of 
boundary layer thickness. In the case of measured diffusivity, the 
effect of turbulent mixing is considered by the assumption that 
thermal diffusivity and kinematic viscosity are specified as 
measured diffusivity. Therefore viscosity and conductivity 
become much larger than molecular diffusivity case. As a result, 
thermal boundary layer thickness becomes much larger and flow 
rate increases dramatically when the diameter increases.   
 
As described above, flow rate is affected by many factors. 
Therefore for the flow rate estimation at all ocean areas, unified 
representation of these factors should be done. Accordingly this 
study introduces new dimensionless number RaR, which was 
proposed by Takahashi et al [3]. The RaR is modification of 
Rayleigh number as shown in eq. (5).  
3
i
R
i
R RgrRa
Rαν
∞ −= ⋅ ,                     (5) 
0
L
i iR dxρ= ∫ ,                        (6) 
0
L
R dxρ∞ ∞= ∫ ,                       (7) 
where ,g α  are gravity acceleration and thermal diffusivity. 
, ,iL ρ ρ∞  are pipe length, inside and outside density. The RaR 
estimates the driving force of this calculation by the density 
difference between inside and outside pipe caused by temperature 
and salinity distributions.  
 
The result of realignment with modified RaR is shown in Figure 9. 
The dimensionless flow rate shown in eq. (8) is used. Aihara et al. 
proposed this number in UWT (Uniform Wall Temperature) 
condition. 
 1G Q
Gr rµ= ⋅ ⋅ ,                      (8) 
where , ,Q Grµ  are flow rate, viscosity and Grashof number. 
Figure 9 also includes the calculation results of Mariana area 
with molecular and measured diffusivity to evaluate the 
possibility of unified representation.  
 
As a result, the relationship between dimensionless flow rate and 
modified RaR has same tendency even in different ocean 
conditions. This means that a unified representation of factors, 
which affect the flow rate, is achieved. As shown in Mariana 
experiment in 2004, the modification of flow rate is necessary in 
small RaR region, however Figure 8 shows that the possibility of 
flow rate estimation in all ocean conditions can be evaluated.  
 
Conclusions 
The numerical simulation for upwelling of deep seawater with 
the perpetual salt fountain is conducted in this study. As a result, 
the velocity profile of the upwelling experiment in Mariana area 
is predicted as M-shape flow and the possibility of reverse flow 
near the axis at middle of the pipe is indicated. Additionally flow 
rate is estimated as 43t/day in the pipe with some modifications. 
Furthermore the upwelling flow in 101.29E and 19.49S ocean 
area near Australia is calculated and flow rate estimation 
depending on ocean conditions is discussed. As a result, it 
becomes clear that the unified representation of ocean conditions 
is achieved by the introduction of new dimensionless number RaR, 
which is modified Rayleigh number, and the dimensionless flow 
rate can be evaluated by modified RaR. 
 
 
   Figure 8. Predicted flow rate in the ocean near Australia 
 
 
 
  Figure 9. Dimensionless flow rate at each condition 
 
References 
[1] H. Stommel, A.B. Arons & D. Blanchard, An 
oceanographical curiosity : the perpetual salt fountain, Deep-
Sea Res., 3, 1956, 152-153. 
[2] K. Tsubaki, S. Maruyama, A. Komiya, H. Mitsugashira, 
Continuous measurement of an artificial upwelling of deep 
seawater induced by the perpetual salt fountain, Deep-Sea 
Res., Part I, 54, 2007, 75-84.  
[3] N. Takahashi, S. Maruyama, S. Sakai, K & Taira, A 
numerical analysis of natural convection using temperature 
and concentration differences, Report of the Institute of 
Fluid Science, 13, 2002, 21-30. 
[4] P. Fofonoff, R. C. Millard Jr., Algorithms for computation of 
fundamental properties of seawater, Unesco technical papers 
in marine science, 1983. 
[5] S. Maruyama, K. Tsubaki, K. Taira & S. Sakai, Artificial 
Upwelling of Deep Seawater Using the Perpetual Salt 
Fountain for Cultivation of Ocean Desert, J. Oceanogr., 60, 
2004, 563-568. 
[6] T. Dylan, Ocean thermal energy conversion: current 
overview and future outlook, Renewable energy, 6, 1995, 
367-373. 
[7] T. Nakasone & S. Akeda, THE APPLICATION OF DEEP 
SEA WATER IN JAPAN, UJNR Technical Report, 28, 1999, 
69-75. 
[8] X. R. Zhang, S. Maruyama, S. Sakai, K. Tsubaki & M. 
Behnia, Flow prediction in upwelling deep seawater – The 
Perpetual Salt Fountain, Deep-Sea Res., Part I, 51, 9, 2004, 
1145-1157. 
397


Numerical Simulation of Upwelling Flow in Pipe Generated by Perpetual Salt Fountain

Abstract

Similar works

Full text

Available Versions

University of Queensland eSpace

CiteSeerX