Evaluation on the tracer simulation under different advection schemes in the ocean model

Advection scheme is one of the core challenges in the computational fluid dynamics, which restricts the capacity of model performance in many research areas including the oceanographic modelling studies. Here in this study, we compared the newly developed algorithm HSIMT (advection scheme with High-order Spatial Interpolation at the Middle Temporal level) with the well-known scheme MPDATA (Multidimensional Positive-Definite Advective Transport Algorithm), in the 1-dimensional idealized simulation and the 3-dimensional realistic river plume simulations. The river plume simulation was done with a mainstream ocean model ROMS (Regional Ocean Modeling System). The results showed that in the 1-dimensional test both HSIMT and MPDATA can converge to the analytical results, but HSIMT converges much faster and does not produce overshooting around the sharp front. Accuracy of HSIMT is also free from the choice of timestep, unlike MPDATA. In the ROMS simulation of a surface-trapped river plume, HSIMT also showed great advantages. Results simulated by MPDATA is highly relied on the model resolution, but when the resolution is high enough the results approached to that simulated by HSIMT. The results of this study could assist the further understanding on capacity of advection schemes and further promote the developments of ocean models.


Introduction
Advection is a key process in ocean dynamics. Transports of salt, temperature, or other contents/properties carried by currents are not only important for themselves, but may also feedback the dynamics through baroclinicity or stratification. Hence, it has been a longterm effort for researchers to introduce physically meaningful, numerically stable, and computationally efficient advection schemes to the models. Unfortunately, few (perhaps none) schemes can precisely simulate the tracer advection process even under very simple and idealized oceanic dynamics, which degrades the fidelity of longterm simulation (Huybers and Wunsch 2010). There is a substantial paradox of advection schemes: accuracy vs. stability, or in another word, dissipation vs. oscillation. Reaching a higher-order (≥2) accuracy is always desired, which however is often at a cost of numerical oscillation. For example, based on the Taylor expansion, the classic Eulerian scheme or Lax-Wendroff scheme should be 2nd-order accurate, yet they have fatal numerical oscillations thus make the model unstable. While the 1st-order upwind scheme is severely dissipative, it is free from numerical oscillation. In the past decades, many works in the field of computational fluid dynamics were attempting to reach the balancing point between accuracy order and numerical stability, for example, the TVD (Total Variation Diminishing) scheme (Harten 1983;Sweby 1984), MPDATA (Multidimensional Positive-Definite Advective Transport Algorithm) scheme (Smolarkiewicz 1983(Smolarkiewicz , 1984, and the third-order, upstream biased scheme (Shchepetkin and McWilliams 1998), among many others. The constructions of these advection schemes are much more complicated than simply doing the Taylor expansion, and the concept of "upwind" is always essential for the solution to be physically meaningful and numerically stable.
Among these, MPDATA is widely used in the mainstream ocean models, such as in ROMS (Regional Ocean Modeling System) and POM (Princeton Ocean Model). In principle, MPDATA is an iterative scheme, which compensates the implicit diffusivity from the 1storder upwind scheme by introducing a pseudo velocity. A great advantage of MPDATA is that it eliminates (but not thoroughly) the numerical oscillation, while at the same time it possesses higher numerical order than the upwind scheme. Since this nature, ocean models associated with MPDATA scheme are often reasonably stable. In spite of its successes, it has also been noted by the users that MPDATA is diffusive when the simulating for a long period. Wu and Zhu (2010) highlighted it by conducting a long-term 1-D simulation that has an analytical solution, and they proposed another advection scheme in that paper, HSIMT (advection scheme with Highorder Spatial Interpolation at the Middle Temporal level). Based on a series of theoretical and realistic simulations, they showed that HSIMT is more accurate, more stable, and more efficient than many other popular advection schemes including the MPDATA. HSIMT has been successfully used and showed great performance in the numerical studies on the Changjiang River plume (Wu et al., Wu and Zhu 2010, Wu et al. 2011, Wu et al. 2014 and now is a standard option in the ROMS.
Here in this study, we focus on the comparison between the HSIMT and the MPDATA, because MPDATA is widely used and has already shown advantages. This study did not attempt to make a comprehensive comparison on all existence advection schemes, because there have been many similar studies. Instead, we focused on the model performance with regard to different spatial resolution and timesteps. In the section2, a brief review of HSIMT was given. Followed in Section 3 was a 1-D numerical test with analytical solution, and a 3-D surface-trapped river plume simulation with a grid-resolution-convergence analysis. Finally, the conclusion was drawn in Section 4.

A review on the HSIMT
The tracer movement in fluid can be described by an advection-diffusion equation. For the purpose of this study, we dropped the diffusion process. For the advection processes, it is useful to begin with the one-dimensional situation that can be described by the following equation: where a is an active or passive tracer, u is velocity. For 1D incompressible flow, u is constant with x, and when the boundary condition was set to zero at −∞ and +∞ the analytical solution for (1) is: where a o (x) is the initial distribution.
(2) means that the tracer distribution remains unchanged except for a uniform shifting in space. This provides an accurate benchmark for validating the performance of different schemes. MPDATA is well documented in extensive literatures (e.g. Smolarkiewicz 1983Smolarkiewicz , 1984, and the readers can refer to that paper. HSIMT is a third-order in space, second-order in time, and finite-volume-based geometric advection scheme (Wu and Zhu 2010). Although its details can be found in that paper, here we provided a more straightforward description. HSIMT was constructed solely based on the integrated form of the advection eq. (1) over [i − 1/2, i + 1/2]Δx and in [n, n + 1]Δt. According to the Newton-Leibniz theorem and the mean value theorem for integral, we have: There are two key points of (3) should be highlighted. Firstly, both a n+1 i and a n i are the mean value in the grid space [i − 1/2, i + 1/2]Δx, instead of the value at the exact position iΔx. When the spatial tracer distribution is not linear, these two are different. In fact, this point is exactly the way advection process is coded in most numerical models, i.e. the local variation of tracer is determined by the total flux into the control cell volume. On the contrary, most algebraic finite difference schemes (such as MPDATA) treated the a n i as the tracer value at iΔx and nΔt, and discrete (1) based (directly or indirectly) on Taylor expansion to achieve certain accuracy orders. Secondly, the cell-interface fluxes u i − 1/2 a i − 1/2 are defined at the time level (n + θ)Δt according to the mean value theorem for integral, not explicitly at nΔt or implicitly at (n + 1)Δt.
So far, no approximation has been made, and (3) describes the original advection eq. (1). For a C-grid ocean model the velocity u is defined at (i − 1/2)Δx, hence, the main task is to determine a n+θ i−1/2 . Assume u keeps constant (or varies linearly thus is taken as 0.5 u n i + u n+1 i ) during [n, n + 1]Δt, according to (2), we have: where λ ≡ u∆t/∆x is the CFL (Courant-Friedrichs-Lewy number) number. Hence, the cell-interface tracer value at (n + θ)Δt can be determined by back-tracking to its original location at nΔt time-level. Assume that the tracer distribution at nΔt is represented by a local quadratic curve g n (to the third-order approximation) over two neighboring if u > 0 for instance. g n has three unknowns that can be determined, since Hence, In practice, θ can be treated as 0.5 for a reasonable approximation. Therefore, the cell-interface flux can be determined to the third-order approximation and a n+1 i can be determined. Since the high-order advection scheme often introduces numerical oscillation, some sort of flux limiter is required in practice to ensure the model stability. A general TVD (total variation diminishing) limiter by Sweby (1984) was found to work very robustly and efficiently. Hence, the final expression of HSIMT scheme is as follows: where As can be seen, the quadratic interpolation is finally represented by the simple coefficient β L or β R , which makes HSIMT to be very efficient in the computation. (4)

Idealized numerical tests
In this part we compare the performance of HSIMT and MPDATA by using both 1D and 3D idealized experiments.

1D advection with incompressible flow
Conducting 1D test to evaluate the advection schemes has been a widely used method since 1) it has an analytical solution (i.e., Eq. (2)) when the fluid is incompressible, 2) it is easy for the readers to repeat the tests themselves. Once an initial tracer distribution was prescribed, this distribution should keep the shape unchanged, according to (2). In many similar studies a constant velocity was set, which restricted the model simulation period since the tracer will move out of the domain when the simulation time is long. It might be too short for model to produce visible difference between numerical schemes. Alternatively, Wu and Zhu (2010) used a spatially uniform periodic current (to mimic the tidal current in the ocean) that enable the 1D tracer simulation to last for an infinite period.
The model domain was 100 km long and the initial passive tracer distribution was set to trapezoidalshaped (thick black line in Fig. 1). Two symmetric fronts were presented with a tracer range of 4.9 and the front width of 5 km. The velocity amplitude is 1 m/s and varies sinusoidally with a period of 12 hr. Similar numerical experiments were done in Wu and Zhu (2010) that included the tests of more numerical schemes, while in this study we presented a more detailed check on the model performance with regard to the grid size and timesteps. Thus, two sets of numerical experiments were conducted under 1) different grid sizes and 2) different timesteps. For the first set, 80 numerical experiments were conducted with the grid size decreased from 2 km to near zero, so as to check the schemes' convergence performance with regard to the grid size. Ideally, the model should give the same result as the analytical solution if sufficiently small grid size was used. In each run of this set, the "CFL" number (calculated with the maximum velocity) was set to 0.5.
For the second set, we fixed the grid size to 500 m (i.e., 1/10 of the frontal width) and varied the CFL number from 0.05 to 1.0. All the experiments were run for 500 periods.

Response to the grid size
With the grid size decreasing, HSIMT converged to the analytical solution very fast, while MPDATA showed notable dissipation (Fig. 1). When the grid size decreased to 500 m, i.e. 1/10 of the front width, the result of HSIMT is already very close to the true value, while MPDATA smoothed the front significantly. Even the resolution approached to 25 m, MPDATA was still unable to capture the sharp front, while HSIMT performed very well. Furthermore, for some grid sizes, MPDATA gave a maximum tracer value even higher than initial maxima.
Numerical accuracy was quantified with the root mean squared error, shown as Fig. 2. The error and grid size were scaled by the initial tracer concentration range and the front width, respectively. Obviously, the error of HSIMT decayed fast with the decrease of grid size, and the use of grid resolution higher than 1/25 of the front scale seems unnecessary. While for MPDATA, the error was much larger and decayed very slowly. To reach the similar accuracy as HSIMT scheme with certain grid resolution, say 1/10 of the front size for example, the MPDATA scheme needs a grid resolution 5 times higher. For 2D case, the grids amount will increase for 25 times. Higher grid resolution requires the use of smaller timestep to ensure the instability. These altogether increase the computational cost significantly.

Response to the timestep
From the view of Taylor expansion, one would expect that the model error decrease with the use of smaller timestep. Whereas, the model results showed an opposite tendency (Fig. 3). For the HSIMT scheme, under the grid size of 500 m, the model error reduced for 0.03 if increased the CFL number from 0.05 to 1.0, while for MPDATA this improvement was 0.07. The reason lies on two-folds. First, although the use of small timestep gives a small truncation error according to Taylor expansion, it requires more timesteps for model to simulate a given period, which increases the accumulative error. Secondly, if the CFL number approaches to one, i.e., a water particle travels exactly the grid size in one timestep, the tracer value a n+1 i equals analytically to a n i−1 for the advection process, thus no accuracy loss occurs. Yet, as the current speed in this 1D test was variable and the CFL was defined under the maximum speed, the error is not zero when "CFL" was one. It is easy to find out that for a constant velocity with unity CFL, both MPDATA and HSIMT, even upwind scheme, give the perfect results. It is therefore recommended that the model should use the timestep as large as possible to increase the accuracy of advection processes.
However, as the timestep in numerical model is also, maybe mainly, determined by many other factors, the CFL number varies strongly both in time and in location for a real simulation, it is desired that the model accuracy is not sensitive to the timestep. From Fig. 3, it is apparent that the accuracy of MPDATA scheme is notably dependent on the timestep. The error increased for 0.07 if the CFL dropped from 1.0 to 0.05, similar to what would happen if doubling the grid size (Fig. 2). On the contrary, the accuracy of HSIMT was less sensitive to the timestep, which is another notable advantage.

Idealized river plume test
We continue to test the advection schemes by simulating a ROMS-default experiment of river plume (the riv-erplume1 option). The domain shape is shown as Fig. 4, that a freshwater (of zero salinity) runoff at a rate of 1500 m 3 s − 1 flows through a 3-km-wide, 7.5-km-long, and 15-m-deep estuary, into a relatively steep shelf with a slope α of 0.0007. The Coriolis parameter was set to 1.0× 10 − 4 , a value at Latitude around 45 o N. No external forcing was superimposed, and the open boundary was set to radiative. The ROMS default setting of open boundary condition for this case, i.e. the periodic condition, is not recommended, since it will cause artificial downshelf ambient barotropic current, which distorts the natural development of a river plume (Matano and Palma 2013). Since there is no analytical solution for this river plume, we conducted the "grid convergence" test similar to the previous 1D test. The ROMS default grid resolution for this case was 1.5 km in cross-shelf direction and 3.0 km in along-shelf direction, respectively. We refined this resolution by a factor of 2, 4, 8, and 16, in both two directions, respectively. The reason doing so is that, by using sufficiently small grid size, the model result should converge to the analytical solution, or in another word the difference equation should approach to the differential equation, if the numerical scheme used is compatible. Compatibility is a basic requirement for any numerical schemes, and it is easy to prove that HSIMT and MPDATA both meet this requirement. This actually has already been observed in the 1-D test. The modeled plumes on 10th day were analyzed. Figure 5 shows the modeled surface plume by using HSIMT advection scheme. A classic plume bulge was formed, with a great portion of freshwater was trapped off the river mouth (Yankovsky and Chapman 1997). The results under different grid size were similar, only with very slight difference. We averaged modeled salinity and velocity every four grids for the higher resolution run, and the result was almost identical to that from a lower resolution simulation. This means that the slight difference among the runs in Fig. 5 was not caused by the numerical truncation error, but simply by the effect similar to the "image definition". To understand this, one could imagine what would happen if projecting an image to a coarse mesh (like digital camera). Apparently, the image will become fuzzy with the minimum increased and maximum decreased. This was consistent to Fig. 2, that HSIMT converged very fast with the decrease of grid size and the truncation error became sufficiently small. Quantitatively, the bulge size was ~ 20 km in the cross-shelf direction, and the normalized grid sizes with this scale were all less than 1/25, which were all sufficiently small as suggested by Fig. 2. While for MPDATA, shown as Fig. 6, the modeled plume showed notable difference under different grid sizes. The plume was fresher and more stratified (not shown) when the resolution was low. This is understandable according to physical nature of river plume and the numerical feature of MPDATA. When the riverine water flows into the ambient sea, its buoyancy is eroded gradually by the vertical mixing, thus the surface salinity increases. While MPDATA is dissipative (relative to HSIMT), the riverine water can migrate numerically across the front without being physically entrained with the high-salinity water beneath, thus artificially remains the low-salinity. The coarser the grid is, the stronger the numerical dissipation is, and thus the lower salinity the plume will be. When the resolution was very high, i.e. the cross-shelf resolution was 93.75 m, visually the plume pattern became similar to that simulated with HSIMT, which was consistent to the 1D theoretical test. Another notable feature was that the MPDATA often gave a zig-zag plume boundary, for example on 10th day when the cross-shelf resolution was 375 m. In fact, this kind of zig-zag existed in simulations with all the grid sizes. The only difference was the timing it occurred.

Freshwater content
To further quantify the model performance, we calculated the freshwater content in each salinity class based on the concept of isohaline coordinate (MacCready et al. 2002). The freshwater volume bounded by isohaline surface s * is calculated as: where F = (s o − s)/s o is the freshwater fraction, s is the salinity, s o is the reference salinity in ambient ocean that (12) V s * = s<s * FdV was 32 in this case. The freshwater content of each salinity class can be calculated as: The results of Fc are show as Fig. 7. For clarity we note each of these experiments as HSIMTi and MPDATA i, where i = 1, 2, 3, 4 signify the grid refinement coefficient 2 i .
It was observed that the HSIMT gave very consistent freshwater content distribution with different resolutions (Fig. 7). The salinity of plume water ranged mainly from 23 psu to 32 psu, while the peak was ~ 26 psu. The HSIMT4 seemed to give a slightly fresher plume. This is because the very high resolution can help to better characterize the low-salinity part inside the bulge (Fig. 5).
(13) Fc(s) = dV (s) ds While for MPDATA the results were completely different. In consistent to Fig. 6, the MPDATA gave a much fresher plume, with the salinity of plume water ranges from ~ 20 psu to 32 psu, a spectrum much wider than that with HSIMT. It was highly grid-dependent, as one can observe significant difference between MPDATA1 and MPDATA4. It seemed that with the resolution increased, the associated freshwater content distribution approached gradually to that with the HSIMT. It can thus be speculated that if the resolution was high enough, the MPDATA would give the same distribution as HSIMT.

Freshwater transport
We also calculated the freshwater transport carried by the plume-induced buoyant coastal current at the position of 100 km. The freshwater transport was normalized by the river discharge, and the resultant percentage is shown as Table 1. HSIMT scheme gave a freshwater transport capability of 53% ~ 58%, if varying the grid size. The transport became larger when the grid size decreased. This is reasonable, because the finer resolution can better resolve the buoyant front, which gave a sharper density gradient in the cross-shelf direction. The buoyant coastal current was driven under the thermal-wind balance; thus, the sharper front could give a larger current.
MPDATA gave a freshwater transport capacity ranging from 82% to 63%, much more uncertain than that with the HSIMT. Basically, the transport capacity decreased and approached to that with HSIMT if refining the grid. As the accuracy of MPDATA is second order (Fig. 2), the finer grid should give a more accurate result. The freshwater transport capacity did not decrease monotonically with grid size, and under MPDATA2 the value dropped to 63%, even lower than that with MPDATA4. Looking at Fig. 6, one can observe that under MPDATA2 the plume front was at a status of numerical oscillation, which introduced error in the simulation of downshelf freshwater transport. Note we have learned in the 1D test that MPDATA is not only dissipative, but also overshooting, and sensitive both to grid size and timestep.

Conclusion
Both 1-dimensional idealized and 3-dimensional realistic numerical simulations showed that HSIMT is more advantageous than MPDATA. In the 1-dimensional test both HSIMT and MPDATA can converge to the analytical results, but HSIMT converges much faster and  does not produce overshooting around the sharp front. Accuracy of HSIMT is also much more independent on the choice of timestep, unlike the MPDATA. These two advantages are resulted from the fact that HSIMT is constructed solely based on physical principles and does not introduced too much mathematical approximations.
In the ROMS simulation of a surface-trapped river plume, HSIMT also showed great advantages. Results simulated by MPDATA is highly relied on the model resolution, but when the resolution is high enough the results approached the that simulated by HSIMT. Because of the good performance of HSIMT in ROMS, it has been used by researchers in the studies of different areas (e.g., Rutherford and Fennel 2018; Jia and Whitney 2019, to name a few). With more and more applications by the community, the advantages and disadvantages of advection schemes can be further revealed, thus can assist the development of more advanced algorithms and might lead to a better solution of this difficult problem in the field of computational fluid dynamics.