Ocean model with adjustable arrangements of discrete variables: application to strong tidal flows and low-salinity water dynamics

The assessment of marine environmental risk necessitates the simulation of a series of phenomena related to the risk as well as a measurement of creatures exposed to the risk. As a practical tool, the simulation is based on the establishment of a numerical ocean model. Although several decades have passed since the numerical model for ocean dynamics has been presented, there remains room for fundamental approaches to refine the method for computing solutions. This paper is a report of the development of a novel algorithm of the model. In this algorithm, discrete variables are positioned in a grid to maximally elicit the advantages of a numerical scheme adopted to each term in the governing equations and simplify the program structure. The implemented program is applied to a tidal flow and riverine buoyant plume in the Hinchinbrook Channel in the eastern coast of the Australian Continent. The computation reproduces the observed strong oscillatory flows and low-salinity water dynamics well. The proposed method is applicable to the movements of pollutants in regions of freshwater influence.


Introduction
Modeling and simulation have gained greater significance in the environmental risk assessment because the ongoing and coming society will demand reliable information on the risk and ensure the ecosystem sustainability and health. The presence of persistent organic pollutants in seawater is one of the risk factors. The ocean model and simulations using it ought to function as practical tools for realizing the environmental risk assessment. The ocean model can serve as this tool after being promoted to efficiently provide solutions that are useful for risk assessment and risk management.
The present study addresses a numerical computation method for simulating the physical aspect, that is, the fluid dynamic aspect, considering the room for improvement in the method for numerically solving a series of oceanic fluid dynamics equations (Navier-Stokes equation, continuity equation, and transport equations of thermohaline fields). The physical module in the ocean model now even requires more updates to be established as practical software.
Coastal oceans, particularly regions near a river mouth, are environmentally the most important regions because they are a pathway of pollutants from the land [1]. They are referred to as the region of freshwater influence (ROFI) [2]. Many previous studies addressed the ROFI from the standpoint of ocean dynamics, such as the growth of the coastal density current driven by the fluid density gradient [3][4][5][6][7][8].
On the contrary, no answer has come against the question of how pollutants move in a ROFI. Clarifying the dynamics of the pollutants and the marine biota exposed to the pollutants in the ROFI is crucial in developing a simulation-based assessment approach of the marine environmental risk. For developing this approach, the eventual form of a simulation model must cover the biological, chemical, and physical aspects of oceanographic phenomena.
To refine the ocean model applicable to ROFI, this paper focuses on how discrete variables in equations (i.e., flow velocity vector variables and other scalar variables) are distributed in a computational grid. The geometric configuration of the grid employed in this study is rectangular, that is, a space is spanned by a Cartesian coordinate system, in which the vertical coordinate of a point is measured as the distance from the still water surface (z-coordinate). The variable distribution can have various patterns in the rectangular grid. Arakawa and Lamb (1977) [9] showed five variable arrangements in a horizontal two-dimensional (2D) grid (hereinafter referred to as the Arakawa arrangements). They placed the components of the fluid velocity vector and the scalar variables on the edges, corners, and/or center in the grid. Since then, almost all ocean models proposed so far selected a single arrangement from the five patterns (e.g., [10,11]).
Aside from the ocean model, more flexible patterns of the variable distribution in a grid were proposed in studies on computational fluid dynamics [12] (this concept is referred to as multi-moment method): a certain variable is placed on both edges and centers to elicit the advantage of a numerical scheme: applications of the constraint interpolation profile-conservative semi-Lagrangian 3 (CIP-CSL3) scheme [13] to a hyperbolic partial differential equation (an advection equation) are its examples (e.g., [14]). Although this method requires more memories to record the calculated values, the enhancement of hardware for computation has made it possible to implement an ocean model code using the multi-moment method.
The idea proposed in the present study is a hybridization of the Arakawa arrangements and the multi-moment method. In the framework of the multi-moment method, the Arakawa grid is categorized into a single-moment method; thus, it would be generally impossible to hybrid these two different methods. However, the aforementioned hybridization becomes realizable by decomposing the governing equations into reduced equations based on the fractional time step approach and by flexibly selecting the Arakawa arrangements that are most suitable to solving each reduced equation, admitting one to build a numerical method that retains the advantages of the selected Arakawa arrangements and the multi-moment method.
This study implements a numerical computation method based on the idea presented above and applies it to a simulation of the flow in and around the Hinchinbrook Channel in the northeastern coast of Australia ( Fig. 1) to assess the performance of the proposed method. Substantially non-stationary phenomena (i.e., very strong flows and large variations in salinity) occur in the channel [15][16][17], which are generally difficult to simulate and, thus, expected to be a good test to see if the proposed method works well. The simulation and observation for this channel were conducted by Tabeta et al. (2002) [18]. The present study presents simulation results using a renewed model.
The remainder of this paper is structured as follows: Sect. 2 describes the proposed numerical computation method; Sect. 3 explains the method application to the Hinchinbrook Channel; Sects. 4 and 5 present the computed results and discuss these through comparisons with in-situ observations; and Sect. 6 draws the conclusions.

Governing equations, boundary conditions, and mode splitting
The present model numerically solves the governing equations for incompressible fluid dynamics on the self-rotating earth. These equations are described in the Cartesian (xyz-) coordinate system, in which x, y, and z denote the eastward, northward, and upward axes, respectively. The vertical axis z measures the distance of a point from the still water surface. The horizontal components of the Navier-Stokes equation are as follows: where, t denotes time; variables (u, v, w) denotes the x-, y-, and z-components of the flow velocity; p denotes pressure; 0 is the reference density of seawater; f is the Coriolis parameter; and K h , K v are the horizontal and vertical viscosity coefficients, respectively.
Assuming that the horizontal scale of the simulation region is much larger than the vertical scale, the vertical component of the Navier-Stokes equation can be written in the following approximated form (hydrostatic approximation): where g is the gravitational acceleration, and = (s, T, p) is the seawater density that depends on the seawater salinity s , temperature T , and pressure. The international equation of state of seawater (EOS-80) was used to determine from s , T , and p.
The equation of ensuring the continuity of the incompressible fluid is expressed as follows: The transport equations of the thermohaline field are where, D h , D v denote the horizontal and vertical diffusion coefficients, respectively. The term "(Source)" in Eq. 5 is the rate of salinity source. This term is modeled as a function of the freshwater input rate. Equation 6 is displayed because it is generally solved; however, in the application of this study, T was assumed to be constant as a reference value, and Eq. 6 was not solved because the salinity field dominantly drives the buoyant river plume dynamics. The boundary conditions for the vertical viscosity were Neumann conditions on both ends. The surface and bottom stresses were computed using the bulk formula. We omitted the wind effects; thus, the wind stresses were assumed to be always zero. The boundary conditions for the horizontal viscosity were Neumann conditions on both ends: no-slip conditions were applied. The kinematic boundary conditions on the bottom and free surfaces are imposed as respectively, where is free surface elevation and the ocean floor's depths are represented as z + H(x, y) = 0 , where H = H(x, y) denotes the heights of the water column with its surface undisturbed. The normal components of the flow velocity on the boundaries between the water and land grids are always set zero (kinematic boundary condition).
The abovementioned primitive equations were split into external and internal modes to separately treat the solutions with fast and slow temporal variations in a manner similar to that in [10]. This technique is useful in assigning suitable numerical schemes to each mode and in efficiently marching the simulation time.

Governing equations of external mode and derivation of hyperbolic-type equations
The primitive forms of the equations governing the dynamics of shallow water on a rotating plane are which are derived by vertically integrating Eqs. 1 and 2 using the kinematic boundary condition on the bottom and free surface and applying the Leibnitz's rule. In (h, U, V) , h is the water column height, U and V are the depth-averaged eastward and northward components of the flow velocity is the height of the sea bottom from the reference level below the sea bottom. E u , E v are the eastward and northward components of acceleration caused by the forces exerted by the bottom and wind stresses and the baroclinic pressure gradients transferred from internal mode calculations.
The primitive form given in Eq. 9 is transformed into a hyperbolic-type equation form by expressing it in the matrix form as follows: Linearly superimposing A and B with weights n x and n y , respectively, a matrix is defined as follows [14]: The 2D problem is split into two one-dimensional (1D) problems, that is, the x-and y-directional problems, which were sequentially solved. Setting n x , n y = (1, 0) gives the hyperbolic-type equation in the x-direction as follows: , n x + n y = 1.
where c ≡ √ gh is the phase speed of the long waves without the effect of earth's self-rotation, and R + x and R − x represent newly defined variables obtained by solving the abovementioned 1D hyperbolic-type equations in the x-direction.
By setting the weights n x , n y = (0, 1) , the following 1D hyperbolic-type equations in the y-direction are derived as follows: where R + y and R − y are newly defined variables.

Numerical scheme
Equations 13 and 14, which are hyperbolic-type ones, can be solved by the semi-Lagrangian scheme to calculate vari- Being regarded as Riemann invariants, these variables were transported on each characteristic curve for each these six variables. The basic concept of the semi-Lagrangian schemes is to identify the transported variables at an upstream point on the characteristic curves [19]. The velocities of the transportations are equal to the eigenvalues of matrices A and , which are the diagonal elements of Δ x and Δ y , respectively. The semi-Lagrangian scheme employed herein is the CIP-CSL3 [13]. When determining one of the variables located at a certain point, this scheme constructs a curve interpolating between that point and the point located the transportation distance for a time step upstream.
In this study, Eqs. 13 and 14 (1D problems) were sequentially solved rather than solving the 2D problem at once. This method for computing solutions enables us to maximally elicit the good property of constructing the interpolation curves using the semi-Lagrangian scheme.
The presence of the external E u , E v and Coriolis forces in the right-hand sides of Eqs. 13 and 14 compels the vari- u to vary while moving on the characteristic curves. These variations can be incorporated into the semi-Lagrangian algorithm by evaluating these variables at the upstream point and approximating the external and Coriolis forces by applying a finite difference scheme. The trapezoidal rule was used in this study.
The discrete form of Eq. 13 is expressed as follows, where the superscripts "*" represents an unknown quantity, and " + ," "-," and "up" represent the upstream quantities transporting in the x-direction with the velocities of u + c , u − c , and u , respectively. Similar to the above, the discrete form of Eq. 14 is written as follows: where " + ," "-," and "up" represent the upstream quantities transporting in the y-direction with the velocities of v + c , v − c , and v , respectively.

Numerical method for solving Navier-Stokes and continuity equations
The governing equations in the internal mode are almost similar to the primitive equations (Eqs. 1-6), except for the temporal evolution equation of the free surface elevation, which is solved by the external mode. Equations 1 and 2 are numerically integrated using the fractional time step approach and the variable arrangements in a grid (Fig. 2).
In this study, two arrangements were selected from the Arakawa arrangements (Fig. 2), namely the collocation arrangement, in which all the variables were positioned at the cell center, and the staggered arrangements, in which the vector variables are positioned on the edge, and the scalar variable is positioned at the cell center. When solving the Navier-Stokes equation, the collocation arrangement was applied to the reduced equations for the advection, Coriolis, horizontal viscosity, and vertical viscosity, while the staggered one was applied to the reduced equations for the pressure gradient. The staggered arrangement was applied for average adjustment (Eq. 31 shown later) and continuity equations.
Each term in these equations (advection, pressure gradient, Coriolis, horizontal viscosity, and vertical viscosity) was sequentially integrated using schemes that suit to solving the reduced equations. The employment of the fractional time step approach permitted the separate implementations of the optimum scheme to efficiently and accurately produce solutions. Table 1 presents the employed schemes.
The reduced equations for the advection terms are which were solved using the CIP-CSL3 scheme to determine the temporary u and v values.
The equations used to solve the reduced equation for the Coriolis terms are which were solved using implicit Runge-Kutta scheme.
Next, the following horizontal 2D diffusion problems (Eqs. 25 and 26) were solved using the forward time and central space finite difference scheme to update u and v: Coefficient K h was determined as a function of the horizontal velocity shears (Smagorinsky parameterization).
This was followed by solving the vertical 1D diffusion problem expressed as and the solutions of which were determined using the implicit Runge-Kutta scheme. Coefficient K v was determined as a function of the Richardson number [20].
The hydrostatic pressures were determined using Eq. 3 and the surface elevations transferred from the external mode calculation, thereby computing temporal values of u and v by solving Eqs. 29 and 30. In this step, the spatial differentiations were evaluated using the central difference scheme to the staggered arrangements of (u, v) on the edges and p at the cell center. The temporal differentiations were approximated by the central difference scheme.  SL  SL  TEC TEC TEC TEC -(y)  SL  SL  TEC TEC SL  SL  TEC TEC -(z)  SL  SL  TEC TEC TEC TEC SL  SL  -Coriolis force  IR  IR  TEC TEC TEC TEC TEC TEC -Pressure gradient  TEC  TEC  CD  TEC TEC CD  TEC TEC -Viscosity  (x and y) FTCS FTCS TEC TEC TEC TEC TEC TEC -(z)  IR  IR  TEC TEC TEC TEC TEC TEC -Average adjustment  TEC  TEC  TZ  TEC TEC TZ  TEC TEC -Continuity equation - The temporary values of (u, v) calculated by the earlier steps and denoted by u ′ , v ′ were further updated using Eq. 31: This treatment plays a role in adjusting vertical averages of (u, v) and ensuring the fluid continuity satisfaction.
The discrete variables (u, v) at different positions (edge or center) in a grid must be always synchronized in a series of the fractional steps. Whenever the variables at one position are updated in a certain fractional step, the variables at the other position are updated in an auxiliary manner. In this study, the time evolution converting (TEC) [12] in which increments of variables that have been already updated are spatially averaged to determine increments of variables that have been left unchanged. When completing (u, v) updates on the edges of all grids, the vertical velocity w on the lower and the upper surfaces in the grids were determined by vertically integrating Eq. 4. The lower limits of the definite integrals were calculated using Eq. 8. The spatial gradients of H(x, y) were computed using the cubic spline interpolation.

Numerical method for solving transport equation for salinity
An approach similar to the one for the internal mode velocity was adopted to solve the transport equation of salinity. The collocation arrangement was used to the reduced equations for the horizontal advection and horizontal and vertical diffusion terms, while the staggered arrangement was applied to the reduced equation for the vertical advection term (Fig. 3). Table 2 summarizes the employed schemes. Equation 5 was split into the following reduced equations: First, Eq. 32 (advection problem) was solved using the CIP method [19], followed by the diffusion problems (Eqs. 33 and 34). Their solutions were determined using the forward time and centered space finite difference scheme and the implicit Runge-Kutta scheme, respectively. Coefficients D h and D v were determined as a function of the horizontal velocity shears (Smagorinsky    [20], respectively. The source term is modeled below. Provided the river flow volumes per unit time denoted by R , the volume of freshwater discharged into a grid with volume ΔV for time step Δt is written as R Δ t . Supposing that the advection and diffusions are omitted, the salinity increment Δ s for Δt is expressed as Under the R Δt∕ ΔV ≪ 1 condition, it is approximated as follows: which gives the following differential equation by taking a limit Δt → 0, This equation was solved using the modified Euler scheme.

Coupling of external and internal modes
The split two modes exchange the variables (Fig. 4) to conform to the theory expressed by the primitive equations. When the time marching of the internal mode was completed, the baroclinic pressure gradients in the internal mode were vertically integrated in every water column, and these results were transferred to E u and E v used in the subsequent external mode. When a set of time marching of the external mode was completed, the surface elevations determined in each step of the external mode were temporally averaged and transferred to the internal mode to compute the pressure gradients.

Topography
Hinchinbrook Channel (146.25°E, 18.50°S) is a narrow channel between the Australian Continent and Hinchinbrook Island (Fig. 1). It is a part of the Great Barrier Reef along the eastern coast of the Australian Continent [15][16][17][18]. It has a water depth of less than 15.0 m around the southeastern exit of the channel, which increases to 200 m at the outer edge of the reef, outside which the East Australian Current flows southward. The Herbert and Seymour rivers flow into the channel near the southeastern exit. A water depth dataset with 250-m horizontal spatial resolution was constructed using a nautical chart (chart No. Aus259, Australian Paper Charts) published by Australian Hydrographic Office (AHO). This resolution was chosen considering the grid fineness for capturing the physics involved in the coastal ocean dynamics and limitation in computing power. A part of Hinchinbrook Channel and the eastern coast off Hinchinbrook Island were included in the simulation domain. Accordingly, the domain was bounded by the eastern (offshore), northern, western (inside the channel), and southern open boundaries: EOB, NOB, WOB, and SOB in Fig. 1.
The simulation domain was discretized vertically by 14 layers. The uppermost layer had a thickness of 0.5 m. As middle layers positioned lower, they were given larger thicknesses. The lowermost layer had a thickness of 5.0 m.
A scheme with no unnatural wave reflectance [14] was applied in the external mode module. Conditions of constant value and zero-gradient were applied for the internal mode module on the four boundaries.
The flow velocity and the salinity were observed with a moored instrument at a site in the channel (hereinafter referred to as Station Green#3, Fig. 1) [18]. The simulated period was 15 days from 00:00 on February 15, 2002.

Tidal forcing
The dynamics of the modeled ocean was driven by incident tidal waves. A time series in the sea surface elevation around Hinchinbrook Island was produced by superposing the eight principal tidal constituents (i.e., M2, S2, N2, K1, a b c d Fig. 4 Schematic for the exchanges of the computed values between the external and internal modes. Δ t e and Δ t i denote the time step of time marching in external and internal mode modules, respectively. a Discrete time evolutions in external mode, b variable transfer from the external to internal mode, c discrete time evolution in the internal mode, and d variable transfer from the internal to external mode K2, O1, P1, and SA), each of which was approximated as sinusoidal functions with amplitudes and phases determined from tidal harmonic constants (AusTides published by AHO). This time series (Fig. 5) was used as the incident tide on the western open boundary.

Modeling of the river discharge
The riverine freshwater discharges were modeled by estimating the source term, as described in Sect. 2.3.2. The flow volume data of Herbert River (Fig. 6) obtained from the Bureau of Meteorology of Australia were employed [ R in Eq. 38]. The flow volumes were assigned to the three discharge points using these data (Fig. 1). Halves of the flow volume of Herbert River were assumed to discharge at the eastern two discharge points. The western discharge point was the river mouth of Seymour River, whose flow volume was assumed to be half that of Herbert River.

Results
Modulated oscillatory flows were produced in the narrow channel (Hinchinbrook Channel; Fig. 7). Its amplitude slowly varied. It gradually lessened for the initial 5.5 days and minimized at Day 5.5. This was followed by gradual growth and a maximization at Day 11.0. The largest flow speed was 2.5 m/s. This variation is responses to the neapspring tidal cycle.
The period of the temporal variation in the flow at St. Green#3 was mostly 0.5 day. The internal velocity magnitude in the surface layer at St. Green#3 was the largest among those in the water column there, while that in the bottom layer was the smallest. The external velocities at St. Green#3 were mostly middle values of the internal ones in the surface and bottom layers.
The numerical computation produced very strong bidirectional flows in the east-west direction inside the channel (Fig. 7). As for the tidal flows around Hinchinbrook Island, Burrage et al. (2002) [16] reported flow velocities at Pith Reef located offshore of Hinchinbrook Island that were measured with an acoustic Doppler current profiler, which revealed an oscillatory variation with a magnitude of nearly 0.2 m/s. In comparison, the amplitude of the nearshore tidal flow is a few times larger. The flow velocities in the east-west direction observed at St. Green#3 (Fig. 7) had a small amplitude (0.5 m/s) and a large amplitude (1.5 m/s) during neap and spring tides, respectively. The  Here we comment on the number of the observation data shown in (Fig. 7), which were insufficient to fully validate the calculated tidal flows owing to several limitations in the measurement work. More frequent data of flow velocity are needed for a better comparison.
Low-salinity water mass began to form inside the channel immediately after the computation set in (Fig. 8). The area of the low-salinity water widened in the east-west direction  Fig. 1) inside the channel. This water eventually went out the channel and transported mainly northward along the coast of Hinchinbrook Island. The salinity at St. Green#3 began to decrease from 33.0 psu at Day 0.4, then became almost zero at Day 4.0 (Fig. 9). Subsequently, the salinity varied between 0.0 psu and 6.0 psu until Day 7.4. At Day 7.4, the salinity exhibited a trend of increase, reached approximately 33.0 psu at Day 7.8 and was then kept a little bit smaller than 33.0 psu with a periodic variation by a few psu.
The salinity around Hinchinbrook Island had been measured by some earlier works (e.g., [16] and [17]) through the conductivity-temperature-depth profiler and remote sensing. The sea surface salinity (SSS) maps created with remote sensing [16] showed effluences of low-salinity water from the proximity of the river mouth of Herbert River to the offshore and to the coast of Hinchinbrook Island, which is north of that river mouth. A similar spatial distribution of the SSS was also presented by Tabeta et al. (2002) [18]. The SSS maps obtained from the computation were qualitatively consistent with the observations. Burrage et al. (2002) [16] reported very large amplitudes of the periodical variation in salinity in Hinchinbrook Channel, which ranged from 10 to 30 psu. In contrast to this, the observed salinity at St. Green#3 in 2002 remained lower than 6.0 psu from Day 0.0 to Day 7.0, then abruptly increased up to 21 psu around Day 7.8. These consecutive lowering and rising variations in the SSS were the result of the increase in the volume of the riverine water amount (Fig. 6) that exceeded the normal level. In the model, the SSS at St. Green#3 began to rapidly decrease at Day 0.4 and became almost zero at Day 4.5 (Fig. 9). On the contrary, it periodically increased to about 5 psu under the tidal flow effect there. These computed results agreed well with the observed ones. The observed increasing trend of the SSS after Day 7.5 was well reproduced by the numerical computation. The surface layer at St. Green#3 was covered again with high-salinity water after the SSS increasing event. The observed SSSs from Day 9.0 to Day 13.5 were within the range of 23-33 psu, which were overall captured by the computation. However, the SSS decrease and increase around Day 10.5 were not simulated by the model. This will be noted later as a future issue.

Discussion
The strong tidal flows inside Hinchinbrook Channel can be interpreted as a 1D shallow water propagation in the east-west direction. The agreement of the calculated flow velocity with the measured one demonstrated the good performance of the external mode module [14] and the appropriate exchanges of the variables between the external and internal modes. Although the internal flow velocities were practically the same as the external flow velocity, by more closely looking into them, the flow speeds in the bottom and surface layers were found to be smaller and larger than the external ones. This can be attributed to two factors: 1) impedance of the bottom friction on the flows above the seafloor; and 2) adjustment of the vertical averages of the internal velocities to the external velocity.
The mode-splitting technique employed in this study required us to handle an algorithm that was more complicated than that without this technique. It required an adequate implementation of the coupling between the external and internal modes. The good agreements of the computation with the observation revealed the right function of the coupling computations.
The almost zero salinity at St. Green#3 was provided by the freshwater discharges from the three river mouths. A few computations with different river flow volumes were performed ( Fig. 9) to estimate their influences because there is uncertainty in the data of the river's flow volume (Fig. 6). Although the small volumes gave a slightly higher salinity, the influences were minor.
The decrease in the river flow volume, which was a recovery to the normal flow volume, was accompanied by the increase in salinity at St. Green#3 to the offshore salinity. The low-salinity water mass formed by the large river flow volume flowed out of the channel and left northward along the Hinchinbrook Island, while the higher salinity water masses alternatively flowed into the channel. In these significant variations in salinity, the horizontal advection terms in Eq. 5 prevailed over the other terms, necessitating an exact computation. The CIP method worked well to accomplish this, albeit the auxiliary treatment needing to be attached to regularly smooth the spatial gradients of salinity and stabilize the computations.
The numerical method proposed herein selected the variable distributions in a grid that enable one to adopt the simplest and most optimal schemes for each term in the Navier-Stokes equations and in the salinity transport, continuity, and average adjustment equations. This idea permits two or more degrees of freedom in selecting the variable distributions in a grid in future ocean models.
The difference between the computation and the observation of salinity at around Day 10.5 needs any model refinement. The observation suggests an intrusion of a low-salinity water mass to St. Green #3; however, the occurrence of this short-term event was not well reproduced by the computation. This problem may be solved by improving the method of making the topographic data and treating the boundary conditions, which will be tackled in the future studies.
Although this study addressed only the physical aspects of the oceanographic dynamics, other aspects, most especially chemical and biological ones, must also be considered when developing a numerical computational approach for assessing the marine environmental risk in a ROFI [2,6,7]. A limitation in the current model is noted here: the hydrostatic approximation applied in this study may be invalid when a finer spatial resolution is demanded for conducting a risk assessment for a practical purpose, to which, a hydrodynamic model may exhibit its feature.

Conclusion
This study developed a numerical computation method for simulating the dynamic phenomenon in the coastal ocean. In our method, the governing equations were split into external and internal modes. The external mode solutions were computed using a semi-Lagrangian scheme, while the internal mode ones were computed using the fractional time step method, in which the governing equations were decomposed into reduced equations. When solving the reduced equations, the discrete variables were adjustably arranged in a rectangular grid, such that the numerical schemes assigned to each of the reduced equations can work well and the entire structure of the algorithm becomes as simple as possible. The implemented program was applied to a simulation of the tidal flows and transport of low-salinity water at Hinchinbrook Channel in Australia. The developed method was validated by comparing the computed flow velocity and salinity with the observed ones. The following results were obtained: • The computed flow velocities agreed well with the observed results, showing that the semi-Lagrangian scheme applied to the external mode went well.
• The collocation and staggered arrangements chosen in the internal mode can appropriately perform in the ocean model. • The calculated results demonstrated that the proposed model could simulate the substantially non-stationary phenomena in the coastal region, that is, the strong tidal flows and abrupt variations in salinity owing to freshwater transport.