Comparison of Physical to Numerical Mixing with Different Tracer Advection Schemes in Estuarine Environments

The numerical simulation of estuarine dynamics requires accurate prediction for the transport of tracers, such as temperature and salinity. During the simulation of these processes, all the numerical models introduce two kinds of tracer mixing: 1) by parameterizing the tracer eddy diffusivity through turbulence models leading to a source of physical mixing and 2) discretization of the tracer advection term that leads to numerical mixing. Physical and numerical mixing both vary with the choice of horizontal advection schemes, grid resolution, and time step. By simulating four idealized cases, this study compares the physical and numerical mixing for three different tracer advection schemes. Idealized domains only involving physical and numerical mixing are used to verify the implementation of mixing terms by equating them to total tracer variance. Among the three horizontal advection schemes, the scheme that causes the least numerical mixing while maintaining a sharp front also results in larger physical mixing. Instantaneous spatial comparison of mixing components shows that physical mixing is dominant in regions of large vertical gradients, while numerical mixing dominates at sharp fronts that contain large horizontal tracer gradients. In the case of estuaries, numerical mixing might locally dominate over physical mixing; however, the amount of volume integrated numerical mixing through the domain compared to integrated physical mixing remains relatively small for this particular modeling system.


Introduction
Estuaries provide critical habitat to sustain marine life. The transport of tracer quantities, such as salinity, sediment, plankton, nutrients, and pollutants, are crucial factors that determine the quality and sustainability of estuaries. The salinity is also an important dynamical variable in estuaries, such that the mixing of saline water with fresh water has important consequences for estuarine circulation and stratification [1]. These interactions of mixing and stratification can depend on several factors such as the size, depth, forcing characteristics, river inflow [2], ocean salinity, tidal forcing and wind forcing [3,4]. This implies that an understanding the processes that cause mixing is a fundamental requirement in the prediction of estuarine dynamics.

Physical Mixing
Physical mixing is introduced through the parameterization of eddy viscosity, typically through turbulence closure models. In RANS based models, the ensemble averaged momentum equation (1) results in a Reynolds stress term that accounts for the fluctuating flow.
where subscript "i" and "j" represent the individual directions, U i and u i are the mean and turbulent component of velocities respectively, ε ijl Ω j is the rotation tensor, P is the pressure, g i is acceleration due to gravity, ρ and ρ 0 are total and reference density respectively, u j u i is the Reynolds stress tensor, and v is the molecular kinematic viscosity. Along with the Reynolds stress term that is introduced through the ensemble averaged momentum equations, the tracer transport equation introduces a turbulent tracer flux term that accounts for fluctuation of tracer quantities in the flow. The transport equation for the tracer quantities is written as: where C and c represent the mean and fluctuating components of any tracer quantity (such as salinity, temperature or suspended sediment), u j c is the turbulent tracer flux, and v θ is the molecular diffusivity. The goal of all the turbulence models is to parameterize the fluctuating terms i.e., the Reynolds stress u j u i and the turbulent tracer fluxes u j c in terms of the mean flow. Establishing such a relationship is referred to as the "closure model". The Reynolds stress term is related to mean flow through eddy viscosity for momentum (K V ) and the turbulent tracer fluxes are related to mean tracer gradients through eddy diffusivity for tracers (K T ). There are different closure models that might depend on the application and computational costs. The most commonly used turbulence closure models involving ocean modeling applications originate from the family of two-equation generic length scale (GLS) models [18]. Different variations in the turbulence closure models use different parameterizations to represent the production and dissipation of turbulence, along with a variety of parameterizations to represent the stability and wall functions. A detailed analysis of the equations involving different turbulence closure models and their performance was conducted by [19][20][21].

Numerical Mixing
The discretization of the tracer transport advection term (second term in Equation (2)) introduces numerical mixing in models. The accuracy of the advection term depends on grid size, choice of advection schemes, and time step size. An increase in grid resolution resolves finer scales; however, the available computational resources often limit the selection of grid size. This issue is further complicated, because an increase in spatial resolution accompanies an increase in temporal resolution to maintain stability [22] in ocean models.
Based on the differences in the formulation of advection schemes, the resulting numerical errors (diffusive and dispersive errors) can vary for a given grid size. Physically, the diffusion errors manifest as smearing of fronts and dispersive errors cause the formation of unrealistic ripples near the fronts. Additional steps, such as the use of limiters, are taken to avoid the formation of spurious wiggles (creation of artificial maximum or minima of tracer) and to maintain nonnegative tracer quantities. At a discrete level, the advection schemes are required to conserve mass. Overall, a high-order accurate, non-oscillatory, and mass conserving advection scheme is desirable for practical applications. More details on the numerical properties that affect advection schemes are available in [23]. The amount of numerical mixing on tracer mixing that is caused by various advection schemes can be estimated through the Taylor series or Von Neumann approach. These analysis methods are difficult to apply when the advection schemes are formulated by using adaptive spatial stencils, limiters, or possess multi-dimensionality properties. Therefore, it is important to find the amount of numerical mixing that is caused by the choice of advection schemes. Other methods, such as the reference potential energy to quantify the spurious mixing in non-isopycnal models, have been utilized [24,25], and found that the increase in horizontal resolution led to a higher fraction of spurious mixing through the contribution of vertical processes. Recently, a discrete variance decay method for quantifying physical and numerical mixing by calculating the loss of kinetic energy could quantify the mixing due to advective and diffusive fluxes [26].
Similar to these aforementioned studies, the goal of this work is to quantify physical and numerical mixing, in a 3-D RANS based ocean model that is applied to an estuarine environment. The estimation of the physical and numerical mixing is based on the work of Burchard and Rennau, 2008. The quantification of mixing allows for understanding the impact of different advection schemes on the tracer distribution in estuarine environments. This article is organized, as follows: the methods describing the numerical model, the choice of different advection schemes and the method to estimate physical and numerical mixing are discussed in Section 2. In Section 3, we use four test cases to quantify physical mixing. In Section 4, we discuss the effects of different advection schemes on mixing from the model and apply the methods to estimate mixing on a realistic application of the Hudson River. Section 5 consists of the summary and conclusions from the current work.

Numerical Model Description
The open-source Coupled Ocean-Atmosphere-Wave-Sediment Transport (COAWST) numerical modeling system [27] is utilized for the current work. The hydrodynamic component of this coupled system is ROMS (Regional Ocean Modeling System), which is a three-dimensional, free surface, terrain-following numerical model that solves finite-difference approximations of the RANS equations while using the hydrostatic and Boussinesq assumptions [28,29] with a split-explicit time stepping algorithm [30,31]. The barotropic and baroclinic momentum equations are solved using different time steps (i.e., a shorter time step for the barotropic component) while assuming the hydrostatic approximation. ROMS employs a two-way, time-averaging procedure for the barotropic mode that satisfies the 3-D continuity equation. It uses a horizontal curvilinear Arakawa C grid and vertical stretched terrain-following coordinates. Momentum, scalar advection, and diffusive processes are represented while using transport equations. The density field is determined from an equation of state that accounts for temperature, salinity, and suspended-sediment concentrations. Further details regarding the numerical algorithms and techniques present in COAWST model are mentioned in [27].

Different Tracer Advection Scheme Choices in COAWST Model
The modeling system allows for choices of many of the model components, including several options for advection schemes of momentum and tracers. We have chosen and used in this work three different tracer advection schemes for comparison that are fundamentally different in their formulation. MPDATA has been derived from Lax Wendroff (LW) family of schemes and it was developed [32,33] by subtracting spatial approximations of higher order error terms from the first order solution, thus obtaining a second order of accuracy. MPDATA is formulated to conserve the sign of the tracers (positive definite) in multi-dimensions introducing cross derivate terms rather than preserving monotone solutions in one spatial dimension.

Third-order Upstream-biased Horizontal Scheme (U3H)
U3H was formulated [34] while using an upstream biased stencil through parabolic interpolation to obtain a third order accurate advection scheme. Unlike MPDATA, the U3H is a directionally split scheme (no cross derivative) terms. A fourth-order centered scheme has been used for vertical advection while the U3H only solves horizontal terms.

High-order Spatial Interpolation at the Middle Temporal Level (HSIMT)
HSIMT was formulated [35] while using a local quadratic distribution of tracer computed at the current time level resulting in a third order accurate advection scheme. The tracer values are constrained to avoid oscillations and preserve monotonicity while using the Sweby's Total Variation Diminishing (TVD) limiter [36]. It is a directionally split scheme (no cross derivative terms) and it maintains the positive definite characteristics of the initial tracer quantity.

Physical Mixing
One of the ways to quantify the mixing in ocean models is to evaluate the tracer eddy diffusivity (K T ) (as described in Section 1) introduced in the tracer transport equation (Equation (2)). However, diffusivity alone may not be a strong indicator of mixing, because a higher value of eddy diffusivity in a homogeneously distributed tracer quantity is not producing increased mixing. Burchard and Rennau 2008 [37] suggest that the decay of mean tracer variance is a far better measure of mixing. When considering a three-dimensional domain of volume V, e.g., the entire estuary, we can decompose the tracer as C = C + C , where angle brackets denote the volume average and prime denotes the deviation from the average, with the volume average of prime vanishing, i.e., C = 0. The equation describing the decay of tracer variance is written as: where the tracer vector and velocity vector can be decomposed as C = C + C , → u = u + u , respectively, and angle brackets denote the volume average while prime denotes the deviation from the average. Additionally, ∇ is the gradient operator, and K is a three-dimensional matrix with K xx , K yy , and K ZZ being eddy diffusivity coefficients in three dimensions and other non-diagonal elements (e.g., K xy ) in K being 0.

of 23
The volume integral is evaluated to remove some of the terms and further simplify (Equation (3)). The vertical integral of the third term on the left hand side divergence of vertical diffusive flux ∂ ∂z K ZZ is zero, recognizing that the vertical turbulent flux K ZZ vanishes at the surface and bottom. Additionally, when Equation (3) is integrated over the volume V, the integral of fourth term on the left side −2C u .∇C amounts to zero. The divergence of horizontal diffusive fluxes contained in the third term on the left side are relatively much smaller than vertical due to the large horizontal scales in estuaries. Similarly, the horizontal components of the right side of Equation (3) are relatively small when compared to vertical mixing based on the aspect ratio in estuaries. With these reductions, Equation (3) simplifies to the form of: By using the divergence theorem, the second term from the LHS in Equation (4) can be converted to represent the advection of tracer variance across the domain (Equation (5)), wheren is the unit normal vector and dS is the area of the domain boundaries.
shows that the change of tracer variance inside a domain is controlled by the net flux of tracer variance across the domain boundaries and the dissipation of variance due to mixing inside the domain. According to [37], the right hand side term provides the mixing in the model, as it represents the action of homogenizing tracer variance through mixing. Thus, the physical mixing from the model is equal to: where D physical is the physical mixing, K T is the vertical eddy diffusivity for the tracer, and ∂C ∂z is the vertical tracer gradient. Note that the physical mixing in Equation (6) only accounts for vertical mixing and it has been approximated after neglecting the horizontal contributions.

Numerical Mixing
Equation (4) provides a relationship between decay of tracer variance and mixing. It should be noted that, theoretically, the mixing term (RHS of Equation (4)) only accounts for prescribed vertical mixing (Equation (5)). However, in numerical models, the discretization process of tracer quantities to solve for the advective terms introduces the other form of mixing i.e., numerical mixing. Therefore, physical and numerical mixing both contribute towards total mixing in numerical models. The amount of numerical mixing can be computed whole using the expression provided by [37] and calculated at each time step by taking the difference between advected square of the tracer and the square of advected tracer. It is equal to: where, D numerical is the numerical mixing, C is any tracer quantity, A{} is the advection operator (depending on the advection scheme), and ∆t is the model time step. One can independently calculate the mixing in each direction and the resulting values can be positive or negative, depending on the advection scheme for a particular grid size. All advection schemes that intend to reduce diffusion by including correction/anti-diffusive terms (Flux Corrector Schemes-FCS schemes add controlled diffusion/anti-diffusion to reduce overshoots) or preserve monotonicity by including limiters (Flux Limited Schemes) cause negative numerical mixing. The volume integral of numerical mixing that was obtained from Equation (7) can be added to the RHS of Equation (4) to obtain the total mixing from the model. For a closed domain, the second term (advection term) of tracer variance on the LHS of Equation (4) goes to zero. Therefore, the volume integral of mean tracer variance decay should amount to the volume integral of physical and numerical mixing.

Model Simulations
The model simulations include four idealized test cases that vary with increasing complexity to compare physical and numerical mixing with different advection schemes. These idealized test cases include: (1) Wind-induced mixed-layer deepening that involves only physical mixing, (2) Horizontal tracer patch propagation that involves only numerical mixing, (3) Lock-exchange representing an estuary without any tidal or river forcing and involves both physical and numerical mixing, and (4) Estuary with tidal and river forcing resulting in both physical and numerical mixing from the model. The last model simulation presented in the section includes the application of a realistic Hudson River estuary.

Test Case 1: Wind-induced Mixed-layer Deepening
This test case is based on wind-induced deepening of a mixed-layer that was described in observations of [38] and analytically by [39], and computationally described by [40,41]. This mixed-layer deepening case provides a simplified analysis of mixing processes in the vertical direction. Because of the mixing in vertical direction and no advection, this test case is used to verify the implementation of physical mixing terms in the model. The tracer in consideration is temperature that is initialized by prescribing a vertical linear variation at the surface to the bottom and starts to mix under the application of surface stress. The tracer only mixes due to the vertical diffusion. The model domain is 500 m long with a width of 100 m and a depth of 50 m. The number of grid spacing in the x, y, and z directions are 10, 5, and 200 respectively, which lead to a vertical resolution of 0.25 m. The north and south boundaries are prescribed with free-slip walls while east and west boundaries are periodic (the longest axis is in the east-west direction). Initially, the temperature varies linearly from 13 • C at the surface to 10 • C at the bottom ( Figure 1a). A surface stress of 0.1 N/m 2 in the x-direction provides the surface friction velocity that creates a periodic east-west flow, leading to mixing in the vertical direction. The simulation is performed at a barotropic time step of 10 s for a period of 30 h. The baseline numerical simulation is performed while using the MPDATA advection scheme using the k − ε turbulence model. The total mixing is computed by the volume and time integral of physical and numerical mixing and comparing with the volume and time integral of mean tracer variance decay (Equation (3), described in Section 2.3). For this test case, physical mixing accounts for all of the mixing and it is equal to the total variance decay ( Table 1). The numerical mixing is zero in this test case because the velocity is spatially uniform, which leads to a spatial gradient equal to zero, and thus there is no contribution of advection terms in this case. Therefore, all of the mixing in the model can be attributed to physical mixing. Identical results are obtained with other two tracer advection schemes (U3H and HSIMT).

Test Case 2: Horizontal Tracer Patch Propagation
The second test case is the horizontal propagation of a tracer patch in one-dimension. Because of no horizontal diffusivity and in the absence of any vertical gradients, the mixing in this test case is entirely numerical. The test case helps in verifying the implementation of numerical mixing terms in the model and in comparing the effect of using different tracer advection schemes in the presence of a sharp front. The test case consists of a 3-D model domain that is 100 km long with a width of 100 m and 0.5 m depth. The resolution in the x-, y-, and z, directions is 500 m, 20 m, and 0.05 m, respectively, with the number of grid spacing in x-, y-, and z-directions as 200, 5, and 10, respectively. This leads to a horizontal along-channel grid resolution of 500 m. The model domain comprises of periodic boundary conditions in north, south, east and west directions and the barotropic time step is 40 s. The simulations are first performed without the presence of the tracer to achieve a steady state horizontal velocity of 0.98 m/s. All the vertical layers advect at the same velocity by the application of a body force term through the water column. After a constant streamwise velocity is obtained, the simulations are restarted with the introduction of a passive tracer patch that varies horizontally between 35 km and 65 km initially.   (Figure 1b). The mixed layer depth is found to be 33 m at the end of 30 h (end of simulation). The efficiency of stratified mixing in estuaries determines the rate at which tracers can be vertically transported. One measure of mixing efficiency is the Flux Richardson (R f ) number [10], which is the ratio of buoyant and turbulence production and it is described as: where, the buoyancy production is defined as B = g ∂ρ ∂z K T ; where g is the gravity constant, ∂ρ ∂z is the vertical gradient of density, and K T is the vertical eddy diffusivity. The production term is defined is the vertical gradient of velocity and K v is the vertical eddy viscosity. Figure 1c shows the Flux Richardson number R f for this test case. It can be noticed that R f is close to zero in the near-surface, weakly stratified layer, and it approaches a constant value slightly less than 0.25 through most of the mixing layer.
The total mixing is computed by the volume and time integral of physical and numerical mixing and comparing with the volume and time integral of mean tracer variance decay (Equation (3), described in Section 2.3). For this test case, physical mixing accounts for all of the mixing and it is equal to the total variance decay ( Table 1). The numerical mixing is zero in this test case because the velocity is spatially uniform, which leads to a spatial gradient equal to zero, and thus there is no contribution of advection terms in this case. Therefore, all of the mixing in the model can be attributed to physical mixing. Identical results are obtained with other two tracer advection schemes (U3H and HSIMT).

Test Case 2: Horizontal Tracer Patch Propagation
The second test case is the horizontal propagation of a tracer patch in one-dimension. Because of no horizontal diffusivity and in the absence of any vertical gradients, the mixing in this test case is entirely numerical. The test case helps in verifying the implementation of numerical mixing terms in the model and in comparing the effect of using different tracer advection schemes in the presence of a sharp front. The test case consists of a 3-D model domain that is 100 km long with a width of 100 m and 0.5 m depth. The resolution in the x-, y-, and z, directions is 500 m, 20 m, and 0.05 m, respectively, with the number of grid spacing in x-, y-, and zdirections as 200, 5, and 10, respectively. This leads to a horizontal along-channel grid resolution of 500 m. The model domain comprises of periodic boundary conditions in north, south, east and west directions and the barotropic time step is 40 s. The simulations are first performed without the presence of the tracer to achieve a steady state horizontal velocity of 0.98 m/s. All the vertical layers advect at the same velocity by the application of a body force term through the water column. After a constant streamwise velocity is obtained, the simulations are restarted with the introduction of a passive tracer patch that varies horizontally between 35 km and 65 km initially. Figure 2 shows the tracer concentration with the three different tracer advection schemes (MPDATA, U3H, and HSIMT) and the exact solution propagation after 9600 s. Out of the three schemes, the U3H shows the closest agreement with the exact solution by maintaining a sharper discontinuity when compared to the other schemes. However, the U3H scheme also shows oscillations around the sharp discontinuities of over-and under-shoots causing negative tracer values. The HSIMT scheme shows a sharper front when compared to the MPDATA. The HSIMT and MPDATA both lack the oscillatory behavior of the U3H, which thereby indicates their monotonicity preserving property. The time and volume integrated values of mixing are compared with the mean tracer variance decay ( Table 1). The results show that the total tracer variance decay is equal to the total mixing (in this case only numerical mixing). The least amount of numerical mixing is obtained from the U3H advection scheme (albeit with non-monotonic behavior at the front), while the MPDATA and HSIMT schemes result in similar numerical mixing.
MPDATA both lack the oscillatory behavior of the U3H, which thereby indicates their monotonicity preserving property. The time and volume integrated values of mixing are compared with the mean tracer variance decay ( Table 1). The results show that the total tracer variance decay is equal to the total mixing (in this case only numerical mixing). The least amount of numerical mixing is obtained from the U3H advection scheme (albeit with non-monotonic behavior at the front), while the MPDATA and HSIMT schemes result in similar numerical mixing. Similar to the definition of vertical eddy diffusivity that is related to the physical mixing in Equation (4), one can compute an effective numerical diffusivity that is based on total numerical mixing and relate it to the horizontal salinity gradient. This helps in quantifying the effect of numerical mixing from different advection schemes where horizontal advection dominates. Similar to the definition of vertical eddy diffusivity that is related to the physical mixing in Equation (4), one can compute an effective numerical diffusivity that is based on total numerical mixing and relate it to the horizontal salinity gradient. This helps in quantifying the effect of numerical mixing from different advection schemes where horizontal advection dominates.
where D numerical is the numerical mixing, K advection-scheme is the effective numerical diffusivity coefficient, and ∂C ∂x is the horizontal tracer gradient which in this case is salinity. At the end of the simulation, the numerical mixing that is obtained from Equation (7) and horizontal salinity gradient integrated over the entire domain are used to get the volume averaged effective numerical diffusivity from Equation (9). It is further divided by a factor of 0.5 u∆x, which results in a non-dimensionalized value (Table 2), where ∆x is the spatial resolution. All the schemes have considerably less diffusion than upwind-differencing. The U3H advection scheme shows the least numerical diffusivity, while the MPDATA and HSIMT schemes show comparable effective numerical diffusivity.

Test Case 3: Lock Exchange Problem
Test case 3 is studied to analyze the mixing properties of numerical schemes with strong horizontal and vertical salinity gradients and shear, but without tidal or river forcing. It is setup by introducing less dense water in the western half of the domain (0 psu) and heavier water (30 psu) on the eastern end reproducing the classic "lock exchange problem". The sharp difference in fluid density leads to a shear flow driven by the internal pressure gradient. This results in the formation of fronts near the surface and at the bottom along with a stratified shear layer. Friction and rotation are set to zero.  (Figure 3b). This identifies the location of the front. To further establish that test case 3 is representative of an estuarine system, the interfacial drag coefficient is computed in the middle of the domain at x = 6.4 km. The interfacial drag coefficient (C Di ) is equal to: where ∆u = (u 1 − u 2 ) 2 and u 1 , u 2 are the velocity in the upper and lower layers of the stratified fluid.
The interfacial drag coefficient in the domain at various depths at the end of the simulation reveals a maximum value of C Di = 1.5 × 10 −4 , which is of similar order to a highly stratified estuarine front [1].
For the baseline case with the MPDATA, the spatial distribution of physical and numerical mixing can be compared in Figure 4a,b. It can be noticed that the physical mixing is higher in the pycnocline (high ds/dz, Figure 4a), while numerical mixing dominates at the fronts (high ds/dx, Figure 4b). The spatial distribution of physical and numerical mixing remains similar with the use of different advection schemes (not shown in the figure). To compare the effects of different advection schemes, numerical mixing varying with depth is observed at the end of the simulation (Figure 5a).
It can be noticed that the U3H shows more negative regions of numerical mixing at the fronts when compared to the other two schemes. This is expected from a higher order scheme, such as the U3H, which is formulated to possess more anti-diffusive properties as compared to the MPDATA and HSIMT schemes. To further understand the numerical properties of the three advection schemes, salinity profiles are obtained from model runs at 1.25 m depth (Figure 5b). The gradient of the front is computed with all three advection schemes. The U3H scheme results in a sharper front when compared to the other two schemes, but it also results in an oscillatory behavior with negative salinity. It can be noticed that the HSIMT shows a sharper front than the MPDATA. In addition, the MPDATA shows a phase lag in capturing the sharp discontinuity as compared to the other two schemes. Table 3 provides a comparison of physical and numerical mixing with different advection schemes. Numerical mixing accounts for majority of total mixing with all the advection schemes. The simulation with the U3H scheme shows the lowest numerical mixing while showing the highest physical mixing when compared to the MPDATA and HSIMT schemes. Based on Equation (9), in test case 3, the U3H advection scheme shows the lowest effective numerical diffusivity, while the MPDATA and HSIMT schemes show comparable effective numerical diffusivity (Table 2).
where ∆ = ( − ) and , are the velocity in the upper and lower layers of the stratified fluid. The interfacial drag coefficient in the domain at various depths at the end of the simulation reveals a maximum value of =1.5 × 10 −4 , which is of similar order to a highly stratified estuarine front [1].
(a) (b) For the baseline case with the MPDATA, the spatial distribution of physical and numerical mixing can be compared in Figure 4a,b. It can be noticed that the physical mixing is higher in the pycnocline (high ds/dz, Figure 4a), while numerical mixing dominates at the fronts (high ds/dx, Figure  4b). The spatial distribution of physical and numerical mixing remains similar with the use of different advection schemes (not shown in the figure). To compare the effects of different advection schemes, numerical mixing varying with depth is observed at the end of the simulation (Figure 5a).  Table 3. Effect of grid resolution on physical and numerical mixing using test case 3, the lock exchange. Comparing physical, numerical mixing with total mixing. Column 1: M phys /M total = Ratio of numerical to total mixing; Column 2: M num /M total = Ratio of physical to total mixing; Column 3: M total = Total mixing (non-dimensionalized by total volume) in in (psu) 2 .    It can be noticed that the U3H shows more negative regions of numerical mixing at the fronts when compared to the other two schemes. This is expected from a higher order scheme, such as the U3H, which is formulated to possess more anti-diffusive properties as compared to the MPDATA and HSIMT schemes. To further understand the numerical properties of the three advection schemes, salinity profiles are obtained from model runs at 1.25 m depth (Figure 5b). The gradient of the front is computed with all three advection schemes. The U3H scheme results in a sharper front when compared to the other two schemes, but it also results in an oscillatory behavior with negative salinity. It can be noticed that the HSIMT shows a sharper front than the MPDATA. In addition, the MPDATA shows a phase lag in capturing the sharp discontinuity as compared to the other two schemes. Table 3 provides a comparison of physical and numerical mixing with different advection schemes. Numerical mixing accounts for majority of total mixing with all the advection schemes. The simulation with the U3H scheme shows the lowest numerical mixing while showing the highest physical mixing when compared to the MPDATA and HSIMT schemes. Based on Equation (9), in test case 3, the U3H advection scheme shows the lowest effective numerical diffusivity, while the MPDATA and HSIMT schemes show comparable effective numerical diffusivity (Table 2). Table 3. Effect of grid resolution on physical and numerical mixing using test case 3, the lock exchange. Comparing physical, numerical mixing with total mixing. Column 1: Mphys/Mtotal = Ratio of numerical to total mixing; Column 2: Mnum/Mtotal = Ratio of physical to total mixing; Column 3: Mtotal = Total mixing (non-dimensionalized by total volume) in in (psu) 2 .

Test Case 4: Estuary Test Case with a Rectangular Geometry
Test case 4 is studied to highlight the impact of salt intrusion with varying tracer advection schemes in an idealized estuary. A similar setup was utilized by [41] to study the effects of using different turbulence closure schemes. The test case involves the simulation of an idealized stratified estuary with a rectangular geometry, along with tidal forcing and river input to provide a time-dependent solution of mixing. The purpose of this test case is to evaluate the effect of mixing from the three advection schemes on salt intrusion length. The effect of earth rotation is not included. The model domain is 100 km long, 0.3 km wide with depths linearly decreasing from 10 m at the ocean (western end) to 5 m at the river end (eastern end). The interior of the domain consists of 200 points in the along-channel direction, seven points in the across channel and 20 points in the vertical. This results in a horizontal resolution of 500 m and a vertical resolution of~0.5 m maximum (varying slightly tidally and decreasing landward).
The tidal velocity on the ocean end is set to sinusoidally vary with an amplitude of 0.4 m/s, while the river inflow velocity at the east boundary is a constant 0.08 m/s. Salinity is forced with a radiation boundary condition and nudged with a three-hour time scale to a value of 30 psu at the ocean end and clamped at 0 psu at the riverine end. Temperature is uniform and constant at 10 • C. The free surface on the ocean end is enforced with Chapman boundary condition and radiation conditions are placed on the baroclinic velocities. The surface elevation at the ocean end is subjected to a tidal forcing of 0.5 m with a 12-hour period. The baseline simulation is performed with the MPDATA advection scheme while using the k − ε turbulence model.   Figure 6a shows the initial salinity distribution with the ocean end on the western boundary with 30 psu and river end on the eastern boundary with 0 psu. A stationary simulation is achieved at the end of 10 days. The mixing estimates are made by utilizing a control volume that extends from x = 15 km to x = 100 km. This avoids the open boundary part of the domain that is most affected during the nudging of salinity at the ocean end. The numerical mixing only constitutes 2% of the total mixing for the MPDATA and HSIMT (refer to Table 1). It is the lowest for the U3H advection (i.e., under 1%). The physical mixing is the highest for the U3H simulation, followed by the HSIMT and the lowest physical mixing is obtained in the simulation with the MPDATA advection scheme. The total mixing from the model (sum of physical and numerical mixing) is the highest for the U3H simulation and lowest for the MPDATA simulation. The impact of mixing on the salt intrusion length during peak flood from different advection schemes can also be noticed (Figure 7). The baseline simulation with the MPDATA shows the least salt intrusion with the isohaline of 2 psu at x = 53 km. The salt intrusion with the HSIMT advection scheme is up until x = 56.5 km, while the use of the U3H advection scheme shows salt intrusion up until x = 58 km. This shows that during peak flood, the MPDATA scheme causes the least amount of total mixing, and it also results in the lowest salt intrusion length. On the other hand, the U3H that is associated with the highest total mixing among the three advection schemes also results in the longest salt intrusion length.

Realistic Application of Hudson River Estuary
The Hudson River is located along the northeast coast of the U.S. The skill assessment of the COAWST model for the 3-D salinity structure of the estuary has been successfully studied in the past [41].

Realistic Application of Hudson River Estuary
The Hudson River is located along the northeast coast of the U.S. The skill assessment of the COAWST model for the 3-D salinity structure of the estuary has been successfully studied in the past [41]. The model grid extends from the south at the Battery to the north in Troy, NY (Figure 8a For the tracer quantities, the lateral boundary conditions consist of closed boundaries on the western edge. The eastern and southern boundaries are radiation nudging, while the northern boundaries are clamped. The input of lateral rivers occurs at several points in the domain. More details on the input from lateral rivers is explained in [43]. The tides are input from the ADCIRC model while using a low pass filter at Kings Point and Sandy Hook. This leads to a tidal range on the order of 2 m and depth-averaged tidal velocities reach 1 ms -1 during spring cycle. We simulated the period from March 25 to July 11, 2005 (111 days). This is done while using the MPDATA advection scheme at a barotropic time step of 6 s using the k − ω turbulence model. A bottom roughness of 0.0025 m is used in the setup and the output is obtained every hour.
The spatial distribution of mixing is assessed by integrating the flood and ebb component over 111 days in the region of interest of the domain (Figure 8b). The ebb constitutes 60% of the total mixing, while the flood consists of 40%. We only focus on the ebb cycle (Figure 9) in the current work, because the ebb cycle consists of a higher amount of mixing; however, both flood and ebb demonstrate similar concepts but different spatial patterns. For the tracer quantities, the lateral boundary conditions consist of closed boundaries on the western edge. The eastern and southern boundaries are radiation nudging, while the northern boundaries are clamped. The input of lateral rivers occurs at several points in the domain. More details on the input from lateral rivers is explained in [43]. The tides are input from the ADCIRC model while using a low pass filter at Kings Point and Sandy Hook. This leads to a tidal range on the order of 2 m and depth-averaged tidal velocities reach 1 ms -1 during spring cycle. We simulated the period from March 25 to July 11, 2005 (111 days). This is done while using the MPDATA advection scheme at a barotropic time step of 6 s using the − turbulence model. A bottom roughness of 0.0025 m is used in the setup and the output is obtained every hour.
The spatial distribution of mixing is assessed by integrating the flood and ebb component over 111 days in the region of interest of the domain (Figure 8b). The ebb constitutes 60% of the total mixing, while the flood consists of 40%. We only focus on the ebb cycle (Figure 9) in the current work, because the ebb cycle consists of a higher amount of mixing; however, both flood and ebb demonstrate similar concepts but different spatial patterns.  Figure 9a shows the averaged salinity distribution during the ebb cycle. The physical mixing integrated over the ebb cycle is the highest in the center of the channel south of the constriction from the George Washington (GW) Bridge, north of the Battery, and a few locations along the edges (Figure 9b). At the constriction from the George Washington (GW) Bridge, the salinity gradient is intensified and fronts are formed as the channel width increases. The increasing of width leads to a steepening of the pycnocline [44], which results in a higher physical mixing in that region. Numerical mixing is highest along the lateral river edges (Figure 9c). The volume integrated ratio of numerical to total mixing is 3.6%. The spatial distribution of mixing can be observed ( Figure 10) in more detail by looking at a lateral cross-section at a location South of the GW bridge near a latitude of 40.85 • for a particular time (02-April at 1700 h) from the ebb cycle. The lateral salinity front corresponds to the region where physical mixing is seen to be the highest in Figure 9b.  Figure 9a shows the averaged salinity distribution during the ebb cycle. The physical mixing integrated over the ebb cycle is the highest in the center of the channel south of the constriction from the George Washington (GW) Bridge, north of the Battery, and a few locations along the edges (Figure  9b). At the constriction from the George Washington (GW) Bridge, the salinity gradient is intensified and fronts are formed as the channel width increases. The increasing of width leads to a steepening of the pycnocline [44], which results in a higher physical mixing in that region. Numerical mixing is highest along the lateral river edges (Figure 9c). The volume integrated ratio of numerical to total mixing is 3.6%. The spatial distribution of mixing can be observed ( Figure 10) in more detail by looking at a lateral cross-section at a location South of the GW bridge near a latitude of 40.85° for a particular time (02-April at 1700 h) from the ebb cycle. The lateral salinity front corresponds to the region where physical mixing is seen to be the highest in Figure 9b. The lateral front shows that high physical mixing (Figure 10a) occurs in the upper portion of the water column. The physical mixing has contributions from both the vertical diffusion and the vertical salinity gradient. The vertical tracer diffusion shows a maximum on the lateral shoal where the stratification is minimal, and the vertical salinity gradient is maximum near the top of the water column. However, the maximum physical mixing is in the region where the combination of terms K T (Figure 10b) and ds/dz maximizes (Figure 10c). The numerical mixing (Figure 10d) maximizes at locations where the horizontal salinity gradient ds/dx maximizes (Figure 10e) i.e., at the lateral fronts. The application of the realistic case of the Hudson River Estuary provides further evidence that is consistent with the simulation of the idealized test cases, which showed that the maximum physical mixing occurred along the density interface, while the maximum numerical mixing was at the lateral edges that are associated with large horizontal tracer gradients. A detailed investigation of the dissipation of salinity variance and its spatial and temporal variation through spring-neap cycles in Hudson River estuary is presented by [45]. Based on the method of total exchange flow method of salt, this study [45] showed that maximum mixing occurs in the presence of intensified turbulence and high salinity gradients, while numerical mixing is small for the Hudson River estuary model. The lateral front shows that high physical mixing (Figure 10a) occurs in the upper portion of the water column. The physical mixing has contributions from both the vertical diffusion and the vertical salinity gradient. The vertical tracer diffusion shows a maximum on the lateral shoal where the stratification is minimal, and the vertical salinity gradient is maximum near the top of the water column. However, the maximum physical mixing is in the region where the combination of terms ( Figure 10b) and ds/dz maximizes (Figure 10c). The numerical mixing (Figure 10d) maximizes at locations where the horizontal salinity gradient ds/dx maximizes (Figure 10e) i.e. at the lateral fronts. The application of the realistic case of the Hudson River Estuary provides further evidence that is consistent with the simulation of the idealized test cases, which showed that the maximum physical mixing occurred along the density interface, while the maximum numerical mixing was at the lateral edges that are associated with large horizontal tracer gradients. A detailed investigation of the dissipation of salinity variance and its spatial and temporal variation through spring-neap cycles in Hudson River estuary is presented by [45]. Based on the method of total exchange flow method of salt, this study [45] showed that maximum mixing occurs in the presence of intensified turbulence

Effect of Grid Resolution on Mixing
In most ocean modeling applications, the users are limited by the choice of grid resolution from the practical constraints of computational costs. While the scope of this study is limited to varying different tracer advection schemes to understand their impact on mixing; in this section, we vary the grid resolution in the horizontal direction for test case 3 to discuss the sensitivity of physical and numerical mixing for three different grid resolutions. The baseline case consists of a horizontal resolution of 50 m, and two more simulations are performed with 25 m and 12.5 m, respectively. Table 3 shows that an increase in grid resolution results in an increase of physical mixing and a decrease of numerical mixing. This is consistent with the idea that a higher grid resolution captures finer scales of turbulence to increase physical mixing, while the increased resolution of salinity fronts in the test case results in lower numerical mixing. The trend of increased physical and decreased numerical mixing with an increased grid resolution holds true for all the advection schemes. A similar sensitivity of mixing components to grid resolution is also seen in the work of Burchard and Rennau, 2008 [37] and Ralston et al. 2017 [46], where a higher resolution grid accompanies higher physical mixing and lower numerical mixing.

Effect of Advection Schemes on Mixing
From the results, it is understood that the spatial distribution of physical and numerical mixing is independent of the choice of advection scheme. The negative numerical mixing acts locally in the simulations and none of the simulations cause a net global negative numerical mixing. The amount of numerical mixing is consistently lowest in all the tests cases while using the U3H advection scheme. Consequently, the U3H captures sharp fronts, which thus results in a better resolution of finer structures that leads to a larger physical mixing. The U3H also generates spurious tracer concentrations with undershooting and overshooting of tracer concentration. These oscillations are a feature (also referred to as Gibbs phenomena, [47]) of non-monotonic schemes such as the U3H. The U3H scheme also has significant occurrence of negative mixing in regions of sharp gradients, as indicated in the lock-exchange example. This indicates that there is a penalty that accompanies the smaller amount of numerical mixing with the U3H. On the other hand, the HSIMT and the MPDATA schemes show a similar amount of numerical mixing in all the test cases. Both these schemes preserve the positive definite characteristics of the tracer. In addition, the use of the MPDATA scheme also causes lower total mixing (physical + numerical) for the idealized two-dimensional (2-D) estuary (Test case 4) and it affects the amount of salt intrusion. With a lower total mixing, the amount of salt intrusion is the least during peak flood when compared with the other two advection schemes. Based on these results, we utilized only the MPDATA scheme to compare the amount and spatial distribution of physical and numerical mixing for a realistic 3-D Hudson River estuary case.
These results indicate that the relative contribution of numerical mixing tends to be small in stratified regimes with active vertical mixing. Using the same analysis method to quantify physical and numerical mixing, the unstructured grid used by [46] predicted that numerical mixing contributed more than half of the total mixing in the highly stratified Connecticut River estuary. The current work analyzes a different physical system of Hudson River estuary and it shows that the numerical mixing for a structured model (i.e., ROMS) with the MPDATA advection scheme is not a significant contributor towards total mixing. This is the case in the mixed layer deepening experiment, the idealized estuary, and the Hudson River simulation, all of which have active vertical mixing that is driven either by wind stress or tidal motion. In contrast, in regimes where vertical mixing is absent or weak and horizontal tracer gradients are strong, numerical mixing tends to dominate. This is illustrated by the tracer patch propagation experiment and the lock exchange case. The dominance of numerical mixing in the lock exchange case is surprising in light of the strong density gradient in the pycnocline, which may be expected to produce significant physical mixing. However, in the absence of external forcing (such as tides or winds), the lock exchange produces a relatively small amount of vertical mixing, and the numerical mixing that is associated with the fronts dominates. Well mixed estuaries would also be expected to have a relatively high contribution of numerical mixing, because of the absence of vertical mixing in the well-mixed regime.
It remains to be determined to what extent numerical mixing might represent unresolved processes in regions of strong gradients. For example, laboratory experiments reveal strong mixing at gravity current fronts [48], which is consistent with the strong numerical mixing in the lock-exchange experiment. The numerical mixing that occurs along the flanks of the estuary in the realistic Hudson simulation might bear some relation to secondary flow-induced mixing that may occur in the actual estuary at scales that are unresolved by the numerical scheme, although this conjecture needs to be tested with higher resolution models.

Conclusions
The current paper estimates the two sources of tracer mixing: physical mixing originating from the turbulence closure model and numerical mixing originating from discretization methods, with a focus on an estuarine environment using the COAWST modeling system. The implementation of the mixing terms is based on the prior work of Burchard and Rennau, 2008 [37]. The formulation of physical and numerical mixing terms in the model is verified in test cases 1 and 2, respectively, by comparing the time and volume integrated total mixing that amounts to total tracer variance decay. Test cases 2, 3, and 4 are used to find the effect of utilizing three different tracer advection schemes on mixing terms. The three advection schemes used in this work include the MPDATA, U3H, and HSIMT. All these schemes have been differently formulated from each other while providing a higher order of accuracy than the first order upwind scheme. The test cases highlight the mixing properties of the advection schemes. The advection scheme causing the lowest numerical mixing (U3H) results in capturing finer spatial structures that lead to a higher physical mixing. Owing to lower numerical mixing, it preserves the sharp fronts with better accuracy than the other two advection schemes. At the same time, the U3H scheme causes undershoots and overshoots of tracer quantities which may cause negative tracer concentrations. This suggests that there is a penalty associated with the U3H, although the net numerical mixing is lowest. The MPDATA and HSIMT schemes both cause a similar amount of numerical mixing while preventing the occurrence of negative tracer values. The spatial distribution of physical and numerical mixing is independent of the advection schemes and it does not coincide. It is shown that the total amount of mixing that can vary with the choice of advection scheme and can consequently impact the length of salt intrusion during peak flood. Physical mixing dominates in the pycnocline while numerical mixing dominates at the horizontal tracer fronts. It is seen through the application of the 2-D idealized estuary and the 3-D Hudson river estuary that numerical mixing might locally dominate over physical mixing. However, for both these cases with open boundaries, the amount of integrated numerical mixing through the domain as compared to integrated physical mixing remains relatively small.
The model source is available for download at https://code.usgs.gov/coawstmodel/COAWST and data sets discussed herein, and their associated metadata, are hosted and publicly available on and accessible via DOI:10.5066/P90KDWTX (Idealized domains [49] and DOI:10.5066/P95E8LAS (Hudson river estuary [50]).

Author Contributions:
The project was carried under the guidance of J.C.W. and W.R.G.; J.C.W. and T.S.K. developed the model simulations and implemented the method to obtain the mixing in the model. X.L. and W.R.G. helped in the mathematical interpretation of the mixing equations when applied in the model. X.L. and T.S.K. performed the model data analysis. H.W. provided help in the implementation and analysis of the HSIMT advection scheme that was added during the course of this project. T.S.K., X.L., J.C.W. and W.R.G. contributed towards the preparation of the manuscript.