Depth-averaged Tidal Flow Simulation by Stabilized Finite Element Formulation in the Ocean Surrounding Indonesia

The final goal of this study is to predict marine environments in the ocean surrounding Indonesia by developing computational techniques based on computational fluid dynamics (CFD). As a preliminary step towards this goal, we have developed methodologies for simulating tidal flows in this paper. The mathematical model of tidal flows is based on depth-averaged 2D shallow water equations (SWEs). We use a stabilized finite element formulation with the Nitsche-type weak imposition of the slip boundary condition to discretize these equations. We propose new techniques to specify appropriate boundary conditions on curved coastlines and open boundaries. We apply harmonic analysis to compact all computational results. We also introduce a one-way nesting procedure to promote a higher accuracy and reflect larger scale effects in subsequent nested meshes. We have applied the present methodologies to a tidal flow simulation in the coastal area around the northern coast of Bali, Indonesia. In this area, we have been observing environmental factors and tidal current velocity in real time from September 2018. We compared the simulation results with the observation results, and discussed the prediction accuracies of the present methodologies.


Introduction
In Indonesia, marine culture and capture fisheries have been growing and creating new employment opportunities in fishery-related industries.However, the country still has a problem with unstable productivity, because there are many environmental problems, such as the rising sea level, the drastic change in weather, and land-based marine pollution.It is imperative to find sustainable uses of marine resources in Indonesia.Maintaining Indonesia's marine resources will be the first step towards stable income not only for fishers, but also for those engaged in all fishery-related industries.The necessity of this has been seen elsewhere in local communities all over Indonesia.To understand marine environments, it becomes increasingly practical to monitor environment factors digitally.Local systems in a few scattered locations have attempted to address the issue, but there is a need for a system that is versatile, effective, and easy to implement.Therefore, our research project has been developing new digital systems to monitor and predict marine environments in Indonesia from 2016 (see https://www.jst.go.jp/global/english/kadai/h2810_indonesia.html).Part of our research theme is to predict behaviors of fluid phenomena by means of computational fluid dynamics (CFD) and develop methodologies for predicting marine environments in the ocean surrounding Indonesia.As a preliminary step towards this goal, we have developed methodologies for simulating tidal flows in this paper.
The mathematical model of tidal flows is based on depth-averaged 2D shallow water equations (SWEs).To discretize these equations, we use a stabilized finite element formulation (1)(2)(3) with the Nitsche-type weak imposition of the slip boundary condition on curved coastlines, which was proposed by Urquiza et al. (4) for Stokes flows.For the specification of the tidal height on open boundaries, we use the NAO.99b global ocean tide model. (5)In addition, for the specification of the tidal current velocity on open boundaries, we use the solutions obtained using the momentum equation of the linearized SWEs to input the gradient of tidal height obtained by NAO.99b.We apply harmonic analysis to compact all computational results.We also introduce a one-way nesting procedure to promote a higher accuracy and reflect larger scale effects in subsequent nested meshes.We have applied the present methodologies to a tidal flow simulation in the coastal area around the northern coast of Bali, Indonesia.In this area, we have been observing environmental factors and tidal current velocity in real time from September 2018.We compared the simulation results with the observation results, and discussed the prediction accuracies of the present methodologies.

Governing equations
We consider that tidal flows are governed by the depth-averaged 2D SWEs derived by integrating 3D incompressible Navier-Stokes equations over the depth while assuming the hydrostatic pressure distribution and neglecting the vertical acceleration.
The depth-averaged continuity and momentum equations of the SWEs can be written as ( ) and the basic configuration and each notation of the SWEs are shown in Fig. 1, where t is time, x = (x 1 ,x 2 ) T is the horizontal Cartesian coordinate, z is the bottom height from the reference height, H is the total water height, ζ = H + z is the water elevation from the reference height, u = (u 1 ,u 2 ) T is the depth-averaged horizontal velocity, g (= 0.980665 m/s 2 ) is the gravitational acceleration, v is the horizontal eddy viscosity, and the others are written as where ε ij is the strain rate, S ci is the Coriolis term, and S si is the surface and bottom shear stress term.Furthermore, f c and f b are the Coriolis and bottom friction parameters, respectively; here, the effect of the surface shear stress is assumed to be zero.
We adopt the Smagorinsky model (6) for the horizontal eddy viscosity v given by ( ) where C s is the Smagorinsky constant and Δ g is the grid length.
The Coriolis parameter f c can be written as where ω (= 7.2921 × 10 −5 rad/s) is Earth's rotation rate and ϕ is the latitude.
We consider a mixed model for the bottom friction parameter f b given by ( ) In that paper, the quasi-linear advective form of the compressible equations is adopted.Therefore, we rewrite Eqs. ( 1) and ( 2) to the quasi-linear advective form given by where Ω is the domain, U is the vector of conservation variables, which is given by and A i , K ij , and R are respectively the advection matrix, the viscous matrix, and the vector of other terms, which are given by The definitions of the domain and boundaries used in this study are shown in Fig. 2. The entire boundary is divided into three parts, which are the Dirichlet, Neumann, and slip boundaries.In Fig. 2, n = (n 1 ,n 2 ) T is the outward unit normal vector to the boundary and t = (n 2 ,−n 1 ) T is the unit tangent vector.
The boundary conditions can be written as on , .

Numerical methods
Let U h and W h be approximate trial and weighting functions based on the linear triangular element.The semidiscrete stabilized finite element formulation of Eq. ( 7) with the Nitsche-type weak imposition of the slip boundary condition can be written as where the superscript h indicates the discrete approximation, τ SUPG is the stabilization parameter proposed by Tezduyar, (3) C pen is the constant penalty parameter, n e1 is the number of elements, n sb is the number of slip boundary segments, and b S h is the length of slip boundary segments.The last four terms of the equation are derived by Nitsche-type weak imposition, which was proposed by Urquiza et al. (4) for Stokes flows.We note that Takase et al. (2) proposed the stabilized space-time finite element formulation of the SWEs including the shock capturing term.In this study, the space-time formulation and shock capturing term are not adopted because target flows can be assumed not to be so extreme.
Equation ( 15) is discretized through integration with time by the generalized-α method, for which the Navier-Stokes equations of incompressible flows was first proposed by Jansen et al. (7) The method is formulated to obtain the second-order accurate approximation in time, and it permits the use of a relatively large time increment without the constraint of the Courant-Friedrichs-Lewy (CFL) condition in the explicit scheme.

Coordinate system
When a target domain is small, the 2D Cartesian coordinate system (CCS) may be useful, but when a target domain is a large ocean, the direct use of the geographic coordinate system (GCS) employing latitude and longitude as coordinate components is more effective.In this study, the shape of Earth is assumed to be a sphere with Earth's mean radius R (= 6.371 × 10 6 m).We adopt a local tangential plane approximation on each element as shown in Fig. 3, where ϕ is the latitude, λ is the longitude, and the position (ϕ 0 , λ 0 ) is the center of each local element.
This local transformation of each coordinate component and its gradient from the GCS to the CCS can be expressed as ( ) ( ) # ,  # Fig. 3. (Color online) Geographic coordinate system and a local plane approximation on each element.
Although the origin of the CCS changes locally, it is inconsequential because direct coordinate components of Eq. ( 16) never appear in the SWEs.

Preprocessing methods
The bathymetries used in this study have been obtained from the General Bathymetric Chart of the Oceans (GEBCO), whose resolution is 30 arc-seconds (about 1 km).The latest version of GEBCO is called the GEBCO_2014 grid. (8)For coastlines, we use the full-resolution Global Self-consistent Hierarchical High-resolution Geography (GSHHG) data. (9)Both GEBCO and GSHHG data on Earth are freely available.
To generate a triangular finite element mesh, we use Triangle, (10,11) which is also freely available on the web (http://www.cs.cmu.edu/~quake/triangle.html).Triangle has an option for refining preexisting meshes with an area constraint imposing each maximum triangle area.In this study, the maximum triangle area e max A at element e of a previous mesh is given by where H 0 is the still water level at the element center, CFL is the target CFL number (CFL = 1), and Δt is the target time increment.The mesh refinement is repeated several times.
To obtain the bottom height, we need to interpolate from the given bathymetric data onto the simulation mesh.The given bathymetric data is GEBCO's 30 arc-seconds grid, which we call the background mesh.When the background mesh is equivalent to or coarser than the simulation mesh, the bilinear interpolation is effective.On the other hand, when the background mesh is finer than the simulation mesh, a more global interpolation technique that would take the whole available information into account to avoid the local effects of bilinear interpolation should be considered.In this study, we adopt the cubic spline kernel function used for the particle method. (12)Figure 4 shows the relationship between the simulation and background meshes.
The bilinear interpolation of the bottom height z A on node A [see Fig. 4(a)] can be written as , , , , ξ ξ ξ = is the local position from x A to x i,j and is the grid length.
The interpolation using the kernel function w(r,h) [see Fig. 4(b)] can be written as ) where h is the parameter of kernel scale (we chose h = 0.5r max , where r max is the length from node A to the closest node), k is the index of the background grid, and n kr is the total number of background grids inside the circle with radius r max .In this study, if n kr ≤ 4, the bilinear interpolation is used; otherwise, the interpolation using the kernel function is used.

Open boundary condition
To obtain the open boundary condition, the NAO.99b global ocean tide model representing 16 major constituents with a spatial resolution of 0.5° is used. (5)In this study, we only use eight major constituents, namely, M 2 , D 2 , N 2 , and K 2 for semidiurnal tides, which are about two cycles per day, and K 1 , O 1 , P 1 , and Q 1 for diurnal tides, which are about one cycle per day.NAO.99b can provide the tidal height that corresponds to the water elevation from the mean water level in this study.Conversely, the tidal current velocity is not provided.Therefore, we have attempted to predict the tidal current velocity on open boundaries using the available NAO.99b data.
The tidal height at the given time t and location x can be calculated as where n tide = 8, ω k is the angular velocity of the major tide k, and the coefficients ( ) ) and b ζk (x) can be obtained from the bilinear interpolation of the gridded data of NAO.99b.The gradient of the tidal height can be calculated as , , where the coefficients c ik (x) and d ik (x) can be obtained from the first-order difference.
The momentum equation of the linearized SWEs can be written as In addition, the tidal current velocity can be written as where the coefficients e ik (x) and f ik (x) are unknown in this stage.The substitution of Eqs. ( 24) and (26) in Eq. (25) leads to the following simultaneous linear equations: Finally, we can obtain the following solutions for the coefficients of the tidal current velocity: , , , .

Harmonic analysis
As the variables of fluids depend on time and location, a substantial number of computational results are obtained.We should consider an extraction technique for the necessary data.In this study, we apply harmonic analysis to compact all computational results.
Let the function f h (x,t) indicate a computational result at the given time t and location x such tidal height and tidal current velocity.The function f h (x,t) can be decomposed into the harmonic function f h (x,t) and its residual v h (x,t) denoted by The harmonic function f h (x,t) can be written as , cos sin , where A 0 (x), ( ) , and b k (x) are the harmonic coefficients.These coefficients can be calculated by minimizing the time-integrated squared residual, which is called the least-squares method.

Nesting procedure
A suitably fine mesh is required in the target coastal area to obtain accurate results.However, the applicability of the present specification techniques to open boundaries is not suitable to such a fine mesh because the grid size of NAO.99b is large.In this study, we apply a one-way nesting procedure in which there are a large-scale domain and also small-scale domains, usually with finer resolution, embedded in the first.The large-scale model runs independently of the small-scale model.The small-scale model, however, extracts information from the large-scale model to provide the open boundary and initial conditions.As all solutions of the large-scale model can be expressed in the harmonic functions as discussed in Sect.2.5, the necessary data can be derived by interpolating the harmonic coefficients from the largescale mesh to the small-scale mesh.

Conditions
A tidal flow in the coastal area around the northern coast of Bali, Indonesia, was simulated.A nested configuration with four levels of nested domains shown in Fig. 5 was implemented.The target domains are the S1 and S2 domains shown in Fig. 5(c).The open boundary condition discussed in Sect.2.5 was only applied to the LL domain.We have two observation stations, namely, ST1 in the domain, located at 8.12805°S, 114.60278°E, and a depth of about 12 m shown in Fig. 5(d), and ST2 in the S2 domain, located at 8.18645°S, 114.82292°E, and a depth of about 20 m.Each current meter was installed at half height of the depth, whereas depths at ST1  and ST2 obtained from GEBCO's bathymetric data were 12.377 and 9.519 m, respectively.The location 7.75°S, 114.75°E shown in Fig. 5(c) is one of the grid points defined at the tidal height of NAO.99b.Table 1 shows computational parameters depending on each domain.Small islands that cannot be resolved were omitted when the finite element mesh was generated.Therefore, the still water depth should be minimally limited to stabilize the solution.To accomplish this, we set the minimum water depth.The minimum element length in Table 1 is the approximate value obtained using the equation The other parameters we set were Smagorinsky constant C s = 0.3, grid length 2 , where A e is the area of the element, bottom friction coefficient C b = 0.0026, and Manning's roughness coefficient n = 0.025 m 1/3 s.We used the free slip boundary condition on curved coastlines and set the constant penalty parameter C pen = 10000, which was already examined using a simple verification problem.The simulation period was 150 d, and the integration period to obtain the harmonic coefficients was 120 d.

Results
We compared the present harmonic analyzed results on the M domain to the results obtained by NAO.99b of eight major constituents with respect to the tidal height at the location 7.75°S, 114.75°E from October 1 to 31, 2018.Figure 6 shows the time history of composed results, and Fig. 7 shows decomposed amplitudes and phases.
The tidal wave tendency of the present result is similar to that obtained by NAO.99b.The amplitudes of diurnal tides such as K 1 and O 1 in the present result are greater than those of NAO.99b.
We compared the present simulation results on the S1 and S2 domains to the observation results with respect to the tidal current velocity at the locations ST1 and ST2 from October 1 to 31, 2018.Figures 8-10 show the time history of current velocity components, the histogram of current velocity magnitude, and a current rose diagram, respectively.The present simulation results only include the tidal current component and its nonlinear dispersion component.On the other hand, the observation results include other components such as wind-driven waves and 3D flows.For this reason, the observation results strongly fluctuate compared with the simulation results.
In Fig. 9, a small current velocity below 1.0 cm/s dominates the present simulation at both ST1 and ST2, and the observation at ST1.For the observation at ST2, the most dominant current velocity ranges from 0.5 to 1.0 cm/s and a high current velocity of more than 6 cm/s is shown, even though it is rare (1.4%).
In Fig. 10, at ST1, the prevailing current direction for the present simulation is SE and that for the observation is WSW or SW, and at ST2, the prevailing current direction for the present simulation is ESE and that for the observation is ENE.This result shows a slight offset between the present simulation and the observation.As mentioned in Sect.3.1, such a slight offset is related to the simulation depth at ST2 being different from the actual water depth.It is considered that the flow in the area very close to the coastline strongly depends on the bathymetric data.If more accurate results are required in the simulation, more accurate bathymetric data should be used.In this study, we used GEBCO's bathymetric data, which only has a spatial resolution of about 1 km.
Figure 11 shows the harmonic analyzed current velocity vector with a background of the bathymetric contour, where Fig. 11(a) is in the S1 domain at the time when the South-North velocity component at ST1 is minimum and the tide is flood, and Fig. 11(b) is in the S2 domain  at the time when the West-East velocity component at ST2 is maximum and the tide is ebb.From Fig. 11(a), it is confirmed that the tendency of the current direction near the coastline is opposite to that of the offshore.It can be found that the tidal current from the offshore may not directly enter the bay.On the other hand, from Fig. 11(b), the current direction is almost eastward, and the bathymetric contour shows that the gradient of the depth around ST2 is very steep.In this area around ST2, coral reefs are the main bottom material; thus, the seabed topography is very complex.Unfortunately, GEBCO's bathymetric data cannot represent it.

Conclusions
In the first half of this paper, we have addressed the methodologies used to simulate 2D depth-averaged tidal flows by a stabilized finite element formulation.Most importantly, we have proposed new techniques to specify appropriate boundary conditions on curved coastlines and open boundaries.By applying harmonic analysis to all computational results, we have reduced the number of computational results that can be expressed using the harmonic coefficients.This technique is also effective for one-way nesting, because the open boundary and initial conditions for the nested mesh can be derived by interpolating the harmonic coefficients from a covered course mesh.Moreover, the harmonic analyzed data can be utilized for an unsteady advection diffusion model, which is the mathematical model of transport and dispersion for environmental factors such as pollutants, nutrients, and tracers.This would lead to the next expansion predicting marine environments.
In the second half of this paper, we have applied the present methodologies to a tidal flow simulation in the coastal area around the northern coast of Bali, Indonesia.For the tidal height, the present simulation results are in rough agreement with the results obtained by NAO.99b; however, diurnal tides are greater than those obtained by NAO.99b.For the tidal current velocity, both simulation and observation results at the observation points are very small.The present simulation results only include the tidal current component and its nonlinear dispersion component.On the other hand, the observation results include other components such as winddriven waves and 3D flows.Particularly, the result of the prevailing current direction shows a slight offset between the present simulation and the observation.Such a slight offset is related to the simulation depth being different from the actual water depth.It is considered that the resolution is insufficient for the bathymetric data used in this paper because of the scale of the target area.Thus, bathymetric surveys are needed in the target area to obtain more accurate results.The future challenge for the field observation is to introduce weather stations to observe wind flows and to examine its influence on wind-driven waves.In addition, the bathymetric survey of the target area should be carried out.Moreover, for the present simulation, the surface shear stress term caused by wind flows should be incorporated into the present method, and the simulation should be carried out considering wind flows and using detailed bathymetric data.

Fig. 1 .
Fig. 1. (Color online) Basic configuration and each notation of the SWEs.

Fig. 6 .
Fig. 6. (Color online) Time history of tidal height at the location 7.75°S, 114.75°E.(above: present result on the M domain; below: result of NAO.99b).
Fig. 10.(Color online) Current rose diagram throughout October 2018, which represents the time-integrated current velocity for each direction, where unit is m.(a) ST1 and (b) ST2.

Table 1
Computational parameters depending on each domain.