Quantification of Physical and Numerical Mixing in a Coastal Ocean Model Using Salinity Variance Budgets

Numerical mixing, the spurious mixing primarily generated by the discretization of advection, is often significant in estuarine and coastal models due to sharp, energetic fronts. We compare on‐ and offline estimates of numerical mixing in a submesoscale‐resolving realistic simulation of the ocean state over the Texas‐Louisiana continental shelf. While offline estimates of numerical mixing differ from online estimates, offline methods may be the only analysis available. We use two methods to estimate numerical mixing offline based on the residuals of the salinity squared s2 and volume‐mean salinity variance s′2 ${s}^{{\prime }^{2}}$ budgets. The s′2 ${s}^{{\prime }^{2}}$ budget overestimates the time‐averaged online numerical mixing by 60% at hourly output. The s2 budget compares poorly due to large truncation errors associated with the tendency and advection terms. The residual of the s2 budget starts to converge to the s′2 ${s}^{{\prime }^{2}}$ budget as output frequency increases to 10 min—an unrealistic frequency for long‐term coastal ocean simulations—but neither method unconditionally converges to the online method and therefore cannot be recommended for generic analysis of numerical mixing. We also investigate the effects of horizontal resolution on numerical mixing using a two‐way nested grid with the online method. The volume‐integrated numerical mixing constitutes 57% of the bulk physical mixing—the mixing prescribed by the turbulence closure scheme—in the coarse model and may exceed the physical mixing by half an order of magnitude. We find numerical mixing is reduced by 35% on average in the nested model, likely due to new dynamical processes that emerge in the nested simulation.

In limited-domain coastal ocean simulations, explicit horizontal mixing is commonly kept as low as possible, ideally with no explicit horizontal mixing, allowing a flux-limiting (e.g., MPDATA, Smolarkiewicz & Margolin, 1998) or upwind (e.g., 3rd order, Shchepetkin & McWilliams, 1998) tracer advection scheme to prevent the accumulation of variance at the grid scale. This has the advantage that mixing, in the form of numerical diffusion, is strongest in the regions where tracer gradients are large but is small elsewhere. Tracer distributions across energetic fronts remain monotonic, and weaker tracer gradients associated with less energetic features remain persistent, resulting in tracer fields that, for example, have the same qualitative structures as remotely sensed sea surface temperature or other optical water properties. As such, advection schemes that can cause significant, but targeted, numerical diffusion are important tools in modern coastal oceanography. However, a disadvantage of this approach is that it is difficult to control and quantify the amount of mixing within the simulation.
Numerical methods have been developed for quantifying numerical mixing online during the model run based on tracer variance dissipation (Burchard & Rennau, 2008;Klingbeil et al., 2014), water mass transformation (Holmes et al., 2021;Lee et al., 2002;Megann, 2018), and energetics (Ilıcak et al., 2012;Petersen et al., 2015;Winters et al., 1995). However, online methods are not always practical because they require modifying a model's source code and rerunning it, or may be unavailable when downloading model output from an external source. Consequently, offline methods for quantifying numerical mixing based on tracer variance budgets have become increasingly popular among estuarine and coastal modelers (X. Li et al., 2018;MacCready et al., 2018;Wang et al., 2017Wang et al., , 2021, who often use salinity as the tracer of interest because it dominates the density structure in many coastal regions. Salinity distributions are also strongly controlled by lateral advection-and thus more susceptible to errors in the advection scheme-since salinity structure in coastal regions is primarily driven by lateral boundary conditions (e.g., rivers) rather than surface boundary conditions (e.g., precipitation and evaporation).
We use a subset of a submesoscale-resolving ocean model of the Texas-Louisiana (TXLA) continental shelf in the northern Gulf of Mexico (nGoM) during summer as a case study, where the density field is dominated by salinity variations due to freshwater input from the Mississippi and Atchafalaya rivers. This region is known to have an energetic eddy field (Hetland, 2017) and intense fronts (Qu, Thomas, Wienkers, et al., 2022). Understanding the circulation and mixing in this region is important for better predictions of along-coast connectivity and constituent transport (Thyng, 2023), seasonal hypoxia (Ruiz-Xomchuk et al., 2021;W. Zhang et al., 2015), and the shelf ecosystem (Fennel et al., 2011). We also want to understand the sensitivity of numerical mixing to resolution, so we use a two-way nested child model with five times the parent resolution to test the sensitivity of numerical mixing to horizontal grid resolution.

Quantifying Physical Mixing
Physical mixing can be defined as the dissipation of salinity variance, χ s , due to molecular processes , which is analogous to the dissipation of turbulent kinetic energy ϵ (MacCready et al., 2018). ϵ represents the rate at which kinetic energy is dissipated by molecular viscosity and has units of velocity squared per unit time (m 2 s −3 ). χ s represents the rate at which salinity variance is dissipated by molecular diffusion (Nash & Moum, 2002) and has units of salinity squared per unit time (g 2 kg −2 s −1 ). The molecular dissipation of salinity variance is assumed to be equal to, on average, the dissipation of Reynolds averaged salinity variance at the scales Journal of Advances in Modeling Earth Systems SCHLICHTING ET AL. 10.1029/2022MS003380 3 of 30 of the turbulent eddies (Burchard & Rennau, 2008;Nash & Moum, 2002;Osborn & Cox, 1972), and may be written as where s is the salinity, u is the 3D velocity vector, an overbar is used to denote a Reynolds average, and a prime is used to denote a perturbation from that average. The physical mixing resolved by a model is quantified as the dissipation of Reynolds averaged salinity variance with the eddy diffusivity κ. Due to the vast differences in the horizontal and vertical scales in the ocean, κ is separated into horizontal (κ H ) and vertical (κ s ) components and thus requires different parameterizations. κ H is parameterized by the horizontal mixing scheme and κ s parameterized by the turbulence closure scheme (Burchard & Rennau, 2008;MacCready et al., 2018).
As discussed above, in many estuarine and coastal models, if κ H is prescribed, it is often done with the smallest possible values to reduce spurious noise in the salinity field. Consequently, the horizontal mixing is often much smaller (Geyer & MacCready, 2014) than the vertical mixing , with some exceptions (e.g., Broatch & MacCready, 2022). While this is different than many larger-scale ocean circulation models, we consider horizontal mixing as part of the physical mixing because it is prescribed, whereas numerical mixing is not something that can be directly controlled in the model setup. Generally, horizontal and vertical down-gradient closure models-including parameterizations of vertical shear mixing, harmonic, and biharmonic lateral mixing-result in a positive definite χ s that represents a sink of salinity variance.

Quantifying Numerical Mixing
Salinity variance budgets provide a means to quantify the bulk physical and numerical mixing within a control volume (X. Li et al., 2018;Lorenz et al., 2021;. Numerical models are designed to conserve salinity, but not salinity variance, therefore numerical mixing is treated as an additional subgrid-scale term and quantified as the residual of the salinity variance budget (MacCready et al., 2018). Salinity variance can be defined with or without reference to a volume mean (Burchard & Rennau, 2008;, with the former denoted hereinafter as the volume-mean salinity variance ′ 2 = ( − ) 2 , where is the volume-averaged salinity 1 ∭ , and the latter denoted as the salinity squared s 2 . For example, consider a three-dimensional control volume with no external sources or sinks of volume-mean salinity variance (i.e., no surface salt fluxes or horizontal diffusive boundary fluxes). The numerical mixing  based on the volume-integrated ′ 2 budget is where  denotes a generic discretization operator, t denotes time, V denotes the domain volume, A l denotes the area of the lateral control volume boundaries, u is the velocity normal to the open boundaries, ̂ is the outward pointing normal vector, and χ s is the physical mixing as defined in the right-hand-side (RHS) of Equation 1. Equation 2 is modified version of Equation 3 of MacCready et al. (2018). We have written the volume-mean salinity variance budget inside of a discretization operator because numerical mixing only appears as a consequence of discretization. As with the physical mixing, the numerical mixing is assumed to be a sink of salinity variance (the reasons for this will become more apparent in Section 2), which is reflected in the sign convention of the RHS terms. Equation 2 demonstrates that in the absence of diffusive fluxes, contributions to numerical mixing arise from the time rate of change of variance within the control volume, the advection of variance through the lateral boundaries, and the dissipation of variance due to turbulent mixing.
A salinity variance budget based on ′ 2 (X. Li et al., 2018) showed that the volume-mean salinity variance may be used as a tracer for the stratification within a control volume, an approach that inspired later studies (Broatch & MacCready, 2022;L. Li et al., 2021;Wang & Geyer, 2018 (Burchard, 2020;Burchard et al., 2019), where s in and s out are the inflowing and outflowing salinities calculated with the total exchange flow (TEF) analysis framework (MacCready, 2011), respectively, and Q r is the river discharge. Under the same assumptions,  showed that M may be quantified directly in terms of the net boundary advection of salinity variance into the estuary. Due to the recent invention of the theory, TEF has only been applied to numerical model output. However, a recent paper by Lemagie et al. (2022) suggested how TEF may be measured with field observations.
Without a robust estimate of for many estuaries and coastal regions, it may be feasible to only estimate the advective boundary fluxes of s 2 . Consequently, there is strong motivation to explore how s 2 and ′ 2 are related. Despite the physical mixing being identical in the s 2 and ′ 2 budgets (since has no spatial gradients), this is not necessarily true for other terms (e.g., advection) in their respective equations. This is particularly relevant for realistic, unsteady flows within control volumes characterized by advective and diffusive lateral boundary fluxes, and surface salt fluxes caused by evaporation and precipitation. Therefore, it is unclear a priori if the s 2 budget may be used to robustly estimate numerical mixing given possible dynamical differences with ′ 2 , especially because they will have different numerical implementations when quantified with a numerical model.
Offline computation of the salinity variance budgets introduces additional discretization errors that are dependent on the model output type (i.e., snapshots or averages), output frequency, and numerical methods used to compute any associated derivatives or integrals. Consequently, offline approximations of the salinity variance budgets are not direct estimates of numerical mixing because they include other sources of error besides the truncation errors associated with the advection scheme, and may be further skewed by other sources of spurious mixing. Nevertheless, we follow Wang et al. (2021) and refer to the offline approximation of the salinity variance budgets for quantifying numerical mixing as the "offline method." Wang et al. (2021) examined the accuracy of the ′ 2 budget in a series of idealized and realistic simulations. They suggested the ′ 2 budget agrees well qualitatively with the online diagnostic method of Burchard and Rennau (2008) when the model output frequency captures typical flow timescales (e.g., a gravity current), but they did not quantify the difference between the on-and offline methods for their realistic simulation.
No studies to date, however, have examined the accuracy of the s 2 budget for quantifying numerical mixing. A heuristic argument can be made that the s 2 budget is subject to larger numerical error than the ′ 2 budget because s 2 is generally much larger than ′ 2 unless the system experiences strong, time-dependent stratification (e.g., a salt-wedge estuary). Therefore, it is useful to investigate whether this error significantly impacts the accuracy of the numerical mixing associated with s 2 and ′ 2 .

Objectives of This Paper
Acknowledging the importance of numerical mixing in estuarine and coastal models, the primary goal of this study is to examine numerical mixing in a domain dominated by internal baroclinic instabilities with an energetic flow field (i.e., (1) Rossby number). Previous studies focused on tidally dominated estuaries (Broatch & MacCready, 2022;MacCready et al., 2018;Ralston et al., 2017;Wang & Geyer, 2018) and coastal ocean (Wang et al., 2021) regimes that were strongly advective, but were also characterized by strong physical mixing. In the nGoM during summer, tides and wind forcing are weak; thus we expect our regime to be strongly advective, but without the associated strong boundary layer mixing that characterized previous studies.
Previous studies (referenced above) have suggested success in calculating numerical mixing using offline methods; accurate offline estimates of numerical mixing would be a powerful tool when online methods are not feasible. We investigate the feasibility of this approach in our nGoM domain by expanding the work of Wang et al. (2021) to assess the accuracy of the volume-integrated s 2 and ′ 2 budgets for quantifying numerical mixing. We use the online method of Burchard and Rennau (2008) as the benchmark to assess the accuracy of the offline budgets because it is dependent on a model's internal timestep and insensitive to model output frequency. Offline methods would be considered robust and convergent if they converge to the online method as model output frequency increases, while maintaining relatively high accuracy at lower output frequency. The s 2 budget has the additional constraint that is must yield similar estimates of numerical mixing as the ′ 2 budget at the same output frequency.

Theory
We derive volume-integrated budget equations for salinity squared s 2 and volume-mean salinity variance ′ 2 . The derivation is similar to Lorenz et al. (2021), who derived s 2 and ′ 2 budgets for a realistic estuarine domain. In the case presented here, the salinity surface and bottom boundary conditions are modified for the model used in this study (Regional Ocean Modeling Systems, ROMS, Shchepetkin & McWilliams, 2005). A rigorous exploration of the differences between the s 2 and ′ 2 budgets in the context of numerical mixing is included.

Salinity Squared Budget
Consider a three-dimensional control volume with four open lateral boundaries and an open boundary at the sea surface. We start with Reynolds-averaged salinity conservation in Cartesian coordinates: where u is redefined as the 3D velocity vector. To derive volume-integrated budgets of s 2 and ′ 2 , we first consider the local salinity boundary conditions. ROMS parameterizes surface volume fluxes due to evaporation and precipitation as surface salt fluxes (i.e., a virtual salt flux), which is to say no water is added or taken away from the domain (Nagy et al., 2020;Roullet & Madec, 2000;Shchepetkin & McWilliams, 2005), so the corresponding surface and bottom boundary conditions are where E and P are the evaporation and precipitation per unit area, respectively, ζ is the free surface elevation, and h is the depth of the ocean bottom below mean sea level. To derive an equation for s 2 and the accompanying boundary conditions, we multiply Equations 3 and 4 by 2s and apply the product rule: where the third and fourth terms on the right-hand-side are the dissipation of salinity variance χ s as defined in Equation 1 (Burchard & Rennau, 2008). To form a volume-integrated budget for s 2 , we integrate Equations 5 and 6 over the domain: where A I is the area of the lateral control volume boundaries, A v is the area of the control volume boundaries at the sea surface, and  2 is the s 2 numerical mixing. As shown in Equation 7, five terms affect the time rate of change of s 2 within a control volume (i.e., tendency), the advective flux of s 2 through the lateral boundaries, the diffusive s 2 flux at the sea surface due to evaporation and precipitation, the horizontal turbulent s 2 diffusion within the control volume, the prescribed physical mixing, and the numerical mixing. Note the s 2 diffusion is expressed as a volume integral instead of a boundary area integral to account for the surface and bottom slopes. The numerical mixing arises when Equation 7 is discretized and is calculated as

Volume-Mean Salinity Variance Budget
Following MacCready et al. (2018) and Lorenz et al. (2021), we define the volume-mean salinity variance where is the volume-averaged salinity 1 ∭ . A local equation for ′ 2 is obtained by subtracting from Equation 3, multiplying by 2s′, and applying the product rule: We derive the surface and bottom boundary conditions by multiplying Equation 4 by 2s′: where we have split the ′ 2 surface boundary condition into contributions from s 2 and similar to Lorenz et al. (2021). Noting that the volume integral of s′ is zero and that has no spatial gradients, volume integrating Equations 9 and 10 yields where  ′ 2 is the ′ 2 numerical mixing. Similar to the s 2 budget, the time rate of change of ′ 2 is controlled by five terms: lateral boundary advection, diffusive fluxes through the vertical boundary, horizontal turbulent diffusion within the control volume, physical mixing, and numerical mixing. As with the s 2 budget,  ′ 2 may be calculated by discretizing Equation 11 and rearranging: Note that the physical mixing χ s is the same as in the s 2 budget, but it is unclear how  ′ 2 relates to  2 . Hereinafter, we drop the  notation used to indicate discretization.

Differences Between s 2 and ′ Budgets
An expression relating  ′ 2 and  2 , making use of the Reynolds decomposition used to define the volume-mean salinity variance, is derived to quantify the differences between the s 2 and ′ 2 budgets. If we redefine the salinity as = + ′ , where the overline and prime retain their previous definitions, a local equation for s in terms of and s′ can be written as

Journal of Advances in Modeling Earth Systems
Multiplying Equation 13 by 2 ( + ′ ) leads to an expanded form of the s 2 equation: We note that any terms containing spatial gradients of vanish from Equation 14. However, we retain them for clarity because ⋅ ∇ 2 does not vanish if the divergence theorem is used to transform the advective fluxes into boundary area integrals. Applying a similar decomposition to the surface and bottom boundary conditions and multiplying by 2 After volume-integrating and rearranging to group like terms together, the s 2 budget in terms of + ′ becomes 2 + 2 ∬ ⋅⏟ ⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ 2 tendency and advection

Cross tendency and advection
In current form, Equation 16 contains several new terms due to the influence of as well as the ′ 2 tendency and advection terms. The physical mixing remains the same for both budgets and the resulting numerical mixing is the same as the residual of the s 2 budget. To quantify the differences in numerical mixing between the s 2 and ′ 2 budgets, we obtain an expression for  2 −  ′ 2 after subtracting Equation 11 from Equation 16 and rearranging: As shown in Equation 17, there are six extra terms that appear when expressing s 2 in terms of and s′: two tendency terms, two advective boundary flux terms, and two turbulent diffusion terms. The first tendency term ( 2 ) describes the evolution of the volume-averaged salinity squared with respect to time, and the second term 2 ∭ ′ is a cross term that represents the volume-averaged salinity multiplied by the tendency of the salinity perturbations. The first advection term 2 ∬ describes the advection of the flow through the lateral boundaries multiplied by the volume-averaged salinity squared, and the second advec tion term 2 ∬ ′ is another cross term that describes the advection of the volume-averaged salinity times the salinity perturbations through the lateral boundaries. Subtracting Equation 11 from Equation 16 will yield a residual that is partially

Journal of Advances in Modeling Earth Systems
SCHLICHTING ET AL. 10.1029/2022MS003380 8 of 30 modulated by the model output frequency for unsteady flows. Additionally, the individual terms in each budget (i.e., s 2 tendency and ′ 2 tendency) may evolve differently depending on how large the extra terms are in Equation 17. For clarity, we note that analytically, the left-hand-side of Equation 17 is equal to zero and arises due to discretization errors during the computation of  2 and  ′ 2 . We examine how these extra terms affect the dynamics of the s 2 budget relative to the ′ 2 budget in Section 4.

Realistic Hydrodynamic Model and Nested Grid
We use the Texas-Louisiana continental shelf model (TXLA), which covers the entire Texas-Louisiana shelf and outer slopes ( Figure 1). The model configuration is similar to  and is reviewed here. The model is a validated, realistic implementation of ROMS configured in this iteration as part of the Coupled-Ocean-Atmosphere-Wave-Sediment Transport modeling system (COAWST ver. 3.7, Warner et al., 2010). The model solves the primitive equations over a curvilinear horizontal grid with terrain-following vertical coordinates (Arakawa & Lamb, 1977;Shchepetkin & McWilliams, 2005;. The model uses 30 vertical layers that are concentrated near the surface and bottom to adequately resolve the boundary layers. The horizontal resolution spans from 0.65 km near the coast to 3.7 km near the outer continental slope, with a mean resolution of 1.57 km in the location of the nested grid. A third-order upwind scheme and a fourth-order centered scheme are used for momentum advection and Multidimensional Positive Definite Advection (MPDATA) is used for tracer advection (Smolarkiewicz & Margolin, 1998). The k − ω scheme is used for vertical mixing (Umlauf et al., 2003;Warner et al., 2005) and horizontal mixing is parameterized along geopotential surfaces with constant horizontal viscosity (5.0 m 2 s −1 ) and diffusivity (1.0 m 2 s −1 ) values, both of which are scaled to the grid size. The model uses a baroclinic timestep of 75 s and a barotropic timestep of 1.875 s. The model provides hourly output and neglects tides because they are weak over the region (DiMarco & Reid, 1998). The model is one-way nested into Global HYCOM Reanalysis to provide open boundary forcing and uses ERA-Interim data sets to provide surface forcing (Dee et al., 2011). Additionally, the model uses streamflow data from nine rivers to provide freshwater forcing: Sabine, San Antonio, Trinity, Brazos, Calcasieu, Lavaca, Nueces, including the Mississippi and Atchafalaya, which provides the bulk of the discharge. The streamflow salinity is set to zero at all rivers and streamflow temperature is estimated using the bulk approach described by Stefan and Preud'homme (1993).
The nested model is configured to be two-way such that variables are exchanged across the open boundary. The bathymetry and vertical grid parameters of the child domain are the same as the parent's because the child domain was created from the parent domain. The bathymetry of the child domain was obtained by linearly interpolating the bathymetry of the parent domain. The ratio of the nesting is five to one, resulting in a mean horizontal resolution of 315 m in the child model. The surface fluxes, open boundaries and initial conditions were also obtained from the parent model by interpolating the parent model outputs to the child grid. By doing so, the models can reduce initialization shocks typically seen at the beginning of the simulation. A comparison between the parent and child grids along the open boundaries (i.e., contact points) showed that the fluxes in and out of the boundaries were internally conserved.
The model has been used in several recent studies focusing on small-scale processes and dynamics in the nGoM (Kobashi & Hetland, 2020;Qu et al., 2021;Qu, Thomas, Wienkers, et al., 2022;Ruiz-Xomchuk et al., 2021). For this study, we focus our analysis to the location of the child grid ( Figure 1). The region is located west of the M/A river discharge points, which contribute to the generation of a baroclinic current over the shelf. The region is often saturated with eddies during summer (Figures 2e and 2f) due to formation of baroclinic instabilities (Hetland, 2017;. A diurnal landsea breeze in near-resonance with the local inertial period forces near-inertial motions (e.g., waves and oscillations) that are perturbed by surface and river forcing (X. Zhang et al., 2009). The volume-integrated flow is strongly time dependent and serves as an excellent test case to investigate the accuracy of the offline methods in a realistic simulation. Previous studies focusing on submesoscale processes in the nGoM Luo et al., 2016) have primarily focused east of the M/A discharge points in the aftermath of the 2010 Deepwater Horizon oil spill. The wave-mean flow interactions in our study region can enhance mixing and lead to the rapid vertical exchange of biogeochemical tracers (e.g., oxygen), which becomes more pronounced as the horizontal resolution increases (Qu, Thomas, Wienkers, et al., 2022). We anticipate significant numerical mixing due to strong salinity gradients (Figures 2a-2d), particularly near the northeastern boundary of the child domain where the Atchafalaya plume is more prevalent.

Simulation Overview
We performed two numerical simulations: one of the native TXLA model without nesting (hereinafter the "coarse" simulation), and the second with nesting turned on (hereinafter the "fine" simulation). Both simulations are analyzed from 3 June to 13 July 2010. Due to file corruption issues during the restart process of the fine simulation, we removed the following times (in UTC) from both the coarse and fine grid simulations when directly comparing the simulations: 17 June 22:30-18 June 19:30, 19 June 14:30-19:30, and 9 July 18:30.
There are several notable mixing events driven by a combination of surface salt fluxes, wind stress, and freshwater input that otherwise contrast with periods of low physical mixing, allowing us to study numerical mixing under a variety of different environmental conditions. Over the inner shelf during the study period, contributions to density variations from variation in the temperature field are small relative to those from salinity variations. We focus our analysis to metrics that influence numerical mixing, in particular different metrics from the velocity gradient tensor normalized by the Coriolis parameter such as vertical relative vorticity ζ/f, horizontal divergence δ/f, horizontal strain rate α/f, and the horizontal salinity gradient magnitude |∇ H s|. They are defined as: To isolate the effects of nesting on model processes, all variables discussed hereinafter from the coarse grid are subsetted offline to the location of the child grid and compared directly with the child grid of the fine simulation. Figure 2 displays a sample of surface fields for the coarse and fine simulations of salinity, |∇ H s|, ζ/f, δ/f, and α/f during a time where numerous eddies are seen over the shelf. These eddies are considered to be submesoscale because they exhibit (1) surface Rossby numbers, as approximated by ζ/f Kobashi & Hetland, 2020;McWilliams, 2016). The salinity of the fine simulation exhibits subtle differences near the northeastern boundary, however the effects of nesting are more striking in the vorticity. Several frontal eddies with anticyclonic cores and cyclonic filaments are resolved by both simulations. The cyclonic filaments are often associated with frontal convergence (Kobashi & Hetland, 2020;Qu, Thomas, Wienkers, et al., 2022) and salinity gradients orders of magnitude stronger than those found in the anticyclonic cores where the salinity is more homogeneous. The horizontal strain rate is enhanced throughout the fine simulation, which acts to sharpen the horizontal salinity gradients and enhance frontogenesis (Hoskins & Bretherton, 1972).  To better understand the impacts of nesting on surface processes, Figure 3 displays probability density functions (PDFs) of the various properties discussed above. Figure 3 also displays the latter three quantities sorted by ζ/f > 1, which correspond to (1) surface Rossby numbers associated with submesoscale fronts. The associated median and median-skewness of the velocity gradient tensor quantities and |∇ H s| are shown in Table 1. The relative vorticity is skewed cyclonically (positively) with a negative median for both simulations, with the median of the fine simulation decreasing by 159% and skewness increasing by 31%. The divergence for both simulations have positive medians (i.e., divergence) and negative skewnesses, with the fine simulation increasing over 260% and the skewness decreasing by 34%. The strain for both models follows a χ distribution, consistent with the results of Shcherbina et al. (2013). The salinity gradient magnitude of the fine simulation has a slightly smaller median and skewness compared to the coarse simulation. When sorted by ζ/f > 1, the latter three quantities have significantly higher probabilities toward the tails of their distributions. This makes intuitive sense because the submesoscale fronts are associated with strong horizontal salinity gradients, elevated convergence/divergence, and elevated strain (McWilliams, 2016). Interestingly, PDFs of |∇ H s| remain essentially unchanged between the coarse and fine simulations, with slightly stronger gradients at the tail of the fine distribution, suggesting that changes to the velocity gradients do not necessarily result in similar changes to the salinity gradients.
A possible explanation for the relatively unchanged salinity gradients in the nested model is due to the horizontal mixing scheme. ROMS scales grid-dependent horizontal diffusivities as where κ 0 is a prescribed background turbulent diffusivity (specified previously) and dA is the lateral grid cell area. The maximum grid cell area is unique to each model grid, resulting in κ H being twice as large on average in the fine simulation compared to the coarse. As we will show in Section 4, this increases the magnitude of the horizontal mixing and turbulent diffusion significantly, which may prevent the horizontal salinity gradients from sharpening as expected.
Trends remain similar for the entire water column (Figure 4), however the distributions of ζ/f and δ/f are slightly more Gaussian. α/f still follows a χ distribution, but is significantly weaker relative to the surface distribution. |∇ H s| changes modestly compared to the velocity gradients, with the coarse simulation having slightly larger salinity gradients near the tails of the distribution.

Implementation of the Online Method for Numerical Mixing
The online numerical mixing  is calculated locally following Burchard and Rennau (2008): where A is the advection operator (i.e., MPDATA) and Δt is the model time step, which yields the numerical mixing in each grid cell. A discretized version of Equation 23 may be found in Burchard and Rennau (2008).
 is computed using s as the representative tracer instead of s′ because s′ requires calculating for the location of the child grid during the model run. An analysis (not shown) of Equation 23 for a first order upwind scheme suggests that  should be identical whether s or s′ is used as the tracer. However, this is not necessarily the case for higher-order, nonlinear advection schemes that employ more sophisticated numerical algorithms. Therefore, it is unclear whether the online method should converge to  ′ 2 in addition to  2 if offline discretization errors are small.
Although we will use the relative agreement between on-and offline methods as the basis to assess this, other sources of spurious mixing such as spurious convection due to high grid cell Reynolds numbers (Ilıcak, 2016;Ilıcak et al., 2012), cabbeling, and numerical diffusion may potentially contaminate both on-and offline methods. Generally, we do not expect these sources to be significant based on the results of Wang et al. (2021). We do not expect spurious convection to be significant because ROMS employs a third-order upstream advection scheme that enforces small grid-scale Reynolds numbers (Ilıcak et al., 2012;Shchepetkin & McWilliams, 2005). Cabbeling is also not likely to be an issue, as work by  slightly east of our study area suggested it is less significant in the M/A plume compared to other oceanic regions. However, we cannot say a priori whether numerical diffusion will be significant. We will revisit these comments in Sections 4 and 5. The numerical implementation of the offline method is shown in Appendix A.

Spatial Structure of the Numerical and Physical Mixing
To motivate our analysis of the physical and numerical mixing, Figure 5 displays cross-sections of  and χ s split into horizontal and vertical components when a strong cross-shelf density gradient is generated by freshwater input for the coarse and fine simulations. A nearshore and offshore pycnocline are observed, where the nearshore pycnocline near 29°N is associated with freshwater input from the Atchafalaya River and the offshore near pycnocline is associated with freshwater input from the Mississippi River (Kobashi & Hetland, 2020). The numerical mixing, although noisy, is approximately 2.75 times larger when averaged over the cross section in the coarse simulation compared to the physical mixing and may be orders of magnitude larger locally. The negative mixing is due to the anti-diffusive properties of MPDATA, which acts to reduce the total mixing but does not have a physical interpretation. The numerical mixing is concentrated primarily where the isopycnals are pinched in the pycnoclines, generally corresponding to strong horizontal salinity gradients. Stronger salinity gradients are located shoreward of 29.25°N, however the numerical mixing is small relative to the offshore pycnocline, likely due to increased grid resolution. Previous studies have shown that numerical mixing is related to the horizontal salinity gradients (Hofmeister et al., 2011;Kalra et al., 2019;Klingbeil et al., 2014;Wang et al., 2021) and we investigate this more in Section 5.
The physical mixing exhibits different spatial trends when broken into horizontal ( ) and vertical components. This is due to two reasons: different parameterizations for the horizontal and vertical turbulent diffusivities and different distributions of the horizontal and vertical salinity gradients. is concentrated at the surface and in the middle of cross-section. The areas of high physical mixing roughly correspond to where the vertical shear and turbulent diffusivity are strong (and large vertical velocity variance). Similar to the numerical mixing, is strongly correlated with |∇ H s| and comprises 9.1% of the total physical mixing ( + ) when averaged over the cross section.
Model nesting produces several notable differences, although the general trends remain the same. There are several locations where the numerical mixing in the fine simulation is stronger than the coarse simulation. The physical mixing averaged over the cross-section in the fine simulation is larger than the coarse simulation, with  now comprising 24% of the total physical mixing. As discussed in Section 3, this is primarily due to an increase in the magnitude of κ H in the fine simulation. The spatially averaged numerical mixing decreases by 51%, with the mean numerical mixing being positive for both simulations. In this particular cross-section, the decrease in average numerical mixing appears to be due to a more symmetric distribution of positive and negative values, which is evident in probability density functions and estimates of skewness over the cross-section. Equation 23 is limited in this application because it does not separate the horizontal and vertical contributions of tracer advection to numerical mixing. Horizontal and vertical tracer advection are often computed with different subroutines in ocean models, which could allow  to be decomposed into components. Such a decomposition could be used to assess whether refinement of the horizontal or vertical grid resolution is required to reduce numerical mixing.
To examine the broader spatial variability, Figure 6 displays depth-and time-integrated χ s and  , the ratio of integrated to , and the ratio of integrated  to χ s . Both and  are strongest near the northeastern boundary due to a large influx of brackish water from the M/A river plume.
is more significant near the southwestern boundary. This is unsurprising because the vertical salinity gradients will be weaker relative to the horizontal further offshore where the plume stratification is weaker (in contrast to the northeastern boundary). The integrated  in the coarse simulation is significantly noisier, with a small patch of negative values concentrated near the northeastern boundary. As a consequence, this will decrease the total mixing and may spuriously alter the timescales of mixing processes associated with submesoscale fronts, although more analysis is needed to confirm this. The integrated  exceeds χ s for a significant portion of the domain in the coarse simulation, with the greatest discrepancy occurring near the southwestern boundary where χ s is weaker.
There is a marked enhancement of integrated χ s in the fine simulation. This is especially true for the ratio of to . As discussed previously, we believe this is due to an increase in κ H of the fine simulation. A newly resolved patch of χ s in the fine simulation spanning almost half the latitudinal extent of the domain is seen east of 93°W.  is reduced throughout the domain of the fine simulation. Additionally, the ratio of  to χ s substantially decreases in the fine simulation, with the numerical mixing exceeding the physical mixing in several patches near the southwestern boundary.

Temporal Variability
Here, we explore the volume-integrated s 2 and ′ 2 budgets and compare the accuracy of the numerical mixing estimated from the ′ 2 budget to the online method. Figure 7 shows a time series for the coarse simulation of surface wind forcing, the terms in Equations 8 and 12, and a comparison between the different methods for quantifying volume-integrated numerical mixing. All wind forcing was spatially averaged over the location of the child grid. The land-sea breeze can generally be seen throughout the simulation, with several exceptions (e.g., 30 June and 7 July) due to transient wind events such as storms.
Time series of volume-averaged s 2 and ′ 2 (Figures 7a and 7b) demonstrate that 2 ≫ ′ 2 and that the budgets exhibit a different time dependency. The volume-averaged ′ 2 rarely exceeds 10 (g kg −1 ) 2 and exhibits reduced inertial variability throughout the simulation. The oscillations displayed in the s 2 and ′ 2 budgets correspond strongly to the local inertial period (∼24 hr) because they are modulated by other variables prone to inertial variability such as the horizontal velocities. Higher frequency oscillations are generated by competition between transient changes in the surface wind field and freshwater input from the M/A rivers.
For s 2 and ′ 2 , the advection and tendency terms are highly correlated and have been combined to reduce clutter. Mathematically, this is expressed using the volume integral of their material derivatives and is written as where c is the tracer. Stronger winds are generally associated with larger tendency and advection terms for s 2 and ′ 2 , which is reflected in their respective material derivatives. The difference between the material derivatives and surface fluxes for s 2 and ′ 2 elucidates the influence that removing the volume average salinity has on the dynamics of ′ 2 (Figures 7d and 7e).
, which we hypothesize is due the larger numerical error associated with the tendency term, and to a lesser extent, the advection term.

Journal of Advances in Modeling Earth Systems
SCHLICHTING ET AL.

10.1029/2022MS003380
16 of 30  Time series from the coarse simulation averaged over the location of the child grid for east-west wind stress τ x , north-south wind stress τ y , and wind speed (a). Volume-averaged s 2 (b), volume-averaged ′ 2 (c), volume-integrated s 2 budget (d) and ′ 2 budget (e), and comparison between the offline numerical mixing and the online numerical mixing (f). The tendency and advection terms of Equations 8 and 12 have been combined using the material derivative ∭ , where c denotes the tracer. The horizontal s 2 and s ′ 2 diffusion terms were omitted from (d) and (e) because they are more than an order of magnitude smaller than the other terms, but are included in their respective numerical mixing calculations.

10.1029/2022MS003380
18 of 30 The tendency term in the s 2 budget is nearly an order of magnitude larger than the ′ 2 budget and will experience larger error because both tendency terms are calculated using the same numerical scheme (i.e., centered finite differences). This error should be reduced by increasing the model output frequency, which we investigate in Section 5. Surface s 2 fluxes exert a much larger influence on the s 2 budget for much of the simulation, with the physical mixing term often remaining the smallest. The opposite is true for the ′ 2 budget, where the term balance often experiences competition between the physical mixing and the material derivative.
Estimates of numerical mixing depend strongly on whether s 2 or ′ 2 is used, with both offline methods demonstrating significant error with respect to time.  2 is significantly noisier than  ′ 2 and may exceed the other methods by over an order of magnitude despite the online method being derived in terms of s and s 2 . Consequently, the s 2 budget should not be used for offline quantification of numerical mixing because  2 does not converge to  ′ 2 for hourly output.
 and  ′ 2 are positive definite, with  ′ 2 almost always overestimating  . The accuracy of  ′ 2 relative to  is strongly time dependent and the time-averaged  ′ 2 is 60% larger than the time-averaged  . From 17-24 June the time-averaged  is 141% larger than the time-averaged  ′ 2 . From 1-8 July however, it is only 12% larger. These results suggest the offline method is not suitable for accurate quantification of numerical mixing in our model. However, it is unclear whether convergence between  2 ,  ′ 2 , and  can be achieved by increasing the output frequency. We discuss possible causes for this, including potential impacts of model output frequency, discretization errors, and other sources of spurious mixing that may contaminate the on-and offline methods in Section 5.
The effects of the extra terms on the behavior of the s 2 budget are further explored in Figure 8. The volume-averaged salinity changes by less than 1 g kg −1 over the simulation and experiences a strong inertial signal until the end of June, after which the inertial signal weakens or vanishes for several days at a time. The differences between the s 2 and ′ 2 budgets are generally dominated by competition between and cross-advective . Physically, this relationship shows that in the presence of stratification, the differences between the volume-integrated s 2 and ′ 2 budgets is partially modulated by the size of the control volume V ( 2 is small relative to V) and the advection of the salinity perturbations through the lateral boundaries. The residual of The horizontal cross diffusion term 2 ∭ ∇ ⋅ ( ∇ ′ ) remains on average several orders of magnitude smaller than the other terms in the s 2 and ′ 2 budgets. This is because the horizontal diffusion terms are small for our model due to prescribed value of κ H . Ultimately, all of the extra terms in the s 2 budget are linked to , with ( 2 ) decreasing in magnitude when amplitude of the inertial variability of decreases. The 2 advection term is not as sensitive to the changes in the inertial variability of and oscillates between first and second order, becoming more important at the end of June as begins to increase.
To investigate the bulk impacts of nesting on mixing, Figure 9 displays time series volume-integrated χ s , , and ratios of fine to coarse χ s ,  , total mixing, and the ratio of χ s and  for both simulations. For the coarse simulation, comprises 2.3% of the total physical mixing in the coarse simulation and 5.8% in the fine simulation. In the coarse simulation, the volume-integrated and time-averaged ratio of numerical to total mixing is 50%. For the fine simulation, this ratio decreases to 32% because the physical mixing increases and the numerical mixing decreases. The time-averaged physical mixing in the fine simulation is 42% larger than the coarse simulation and is almost always larger instantaneously, with the greatest disparity occurring during the wind-driven mixing event near 6 July where the fine simulation physical mixing is a factor of four larger than the coarse simulation. The time-averaged numerical mixing in the fine simulation is 35% smaller on average relative to the coarse simulation and only grows larger than the coarse numerical mixing at several times. The total mixing in the fine simulation is 14% larger on average than the coarse simulation, but exhibits large temporal variability and may be twice as large or reduced by a third compared to the coarse simulation. During the first 3 weeks of June, the total mixing in the fine simulation may be more or less than the coarse simulation due to reduced numerical mixing, but for the remainder of the simulation the relative increase in physical mixing generally exceeds the decrease in numerical mixing. A striking comparison is the ratio of numerical to physical mixing. The numerical mixing is larger than the physical mixing in the coarse simula- 2 tendency and advection (b), 2 ′ tendency and advection (c), turbulent horizontal cross diffusion (d), extra s 2 surface fluxes (e), and differences between the s 2 and ′ 2 numerical mixing  2 −  ′ 2 (f).  2 −  ′ 2 yields a large residual that should decrease as model output frequency increases.

Discussion
We find that numerical mixing constitutes a significant fraction of the total mixing-the sum of physical and numerical mixing-and exceeds the physical mixing for much of the simulation due to strong lateral salinity gradients generated by freshwater input from the Mississippi and Atchafalaya rivers. for the coarse and fine simulations (a). All plots below are derived from these values, displaying fine/coarse physical mixing (b), fine/coarse numerical mixing (c), fine/coarse total mixing (d), and numerical/physical mixing for both models (e). Missing values seen in the time series correspond to output lost during the restart process of the fine simulation and were removed from the parent simulation offline for direct comparison. Note each plot is on a log scale of base 10 or base 2.

10.1029/2022MS003380 21 of 30
Our analysis of the offline tracer budgets suggests the s 2 budget may overestimate the online numerical mixing by over an order of magnitude (Figure 7). Building on the work of Wang et al. (2021), we hypothesized that the differences are primarily due to errors associated with the tendency term of each tracer and can be reduced by increasing the model output frequency. To test our hypothesis and examine whether  2 will converge to  ′ 2 , we performed two additional numerical simulations of the coarse model with output frequencies of 30and 10 min (twice and six times finer than the native resolution, respectively). The results of the simulations are summarized in Figure 10 and discussed in Sections 5.1 and 5.2.
Additionally, the local and volume-integrated numerical mixing constitute significant fractions of the total mixing, even when the model resolution is increased to (300 m) in a two-way nested simulation. Despite the time-averaged spatially integrated numerical mixing decreasing in the fine simulation by 35%, Figure 5 suggests that numerical mixing may be occasionally enhanced at the grid-scale. Numerical mixing is dependent on a number of factors, including the magnitude of the horizontal salinity gradients, grid resolution, model timestep, and the representation of dynamical processes. Wang et al. (2021) suggested that numerical mixing is proportional to the square of the horizontal salinity gradients, however they did not demonstrate that this relationship is robust for higher order and more complex advection schemes. We find, based on Figures 3 and 4 that, on average, the salinity gradients do not increase significantly in the fine simulation, but there are more instances of high values of relative vorticity, divergence, and strain. Even when considering the increased horizontal mixing in the fine simulation, the results suggest that different physical processes have emerged in the fine simulation. We discuss the relationship between horizontal resolution, numerical mixing, and model processes in Section 5.3.

′
As the output frequency increases,  2 beings to converge toward  ′ 2 but still remains noisy, with the two in much better agreement for the 10 min output simulation. However, we find this insufficient because 10 min output is impractical for long-term coastal ocean simulations at this time. To test our hypothesis that the higher inaccuracy associated with  2 is primarily due to errors associated with the tendency term ∭ V ∂ t s 2 dV, we where Δ = γ 30 − γ 10 , with γ denoting the terms in the volume-integrated s 2 budget and the subscripts referring to the model output frequency in minutes. The difference in the estimate of numerical mixing due to differences in model output frequency is given by Δ 2 . We computed the covariance between each term in Equation 25 and Δ 2 divided by the variance of Δ 2 to test the significance of each term, which we denote by the variable q. This provides an estimate of the fraction of the variance in Δ 2 explained by each term. For −ΔTendency, q = 1.271 and for −ΔAdvection, q = −0.270, indicating that ΔTendency will over-predict Δ 2 , whereas −ΔAdvection will compensate for this over-prediction. When combined, q = 1.001, indicating that ΔPhysical Mixing, ΔSurface s 2 flux, and ΔHorz. Diff. do not significantly contribute to Δ 2 , which is because they are orders of magnitude smaller relative to ΔTendency and ΔAdvection. ΔAdvection is also significant, but to a lesser extent than ΔTendency because the horizontal velocities are highly unsteady and therefore sensitive to model output frequency. The physical mixing is computed online and remains more periodic relative to ΔTendency and ΔAdvection, so it is less sensitive to model output frequency. Likewise, the freshwater flux (E − P) has an temporal resolution of several hours that is linearly interpolated prior to the calculation of ΔSurface s 2 flux, so it is unsurprising it remains insignificant. ΔHorz. diff. does not change signficantly because the horizontal s 2 diffusion is well under an order of magnitude smaller than the other terms in the volume-integrated s 2 budget. These results were consistent with scatter plots of Δ 2 and the RHS terms of Equation 25 (not shown).

Convergence Between 
′ and  For 10 min output,  ′ 2 becomes slightly less noisy but the general structure remains unchanged and unconditional convergence (i.e., for all times) is not achieved. Interestingly, convergence does not substantially improve as output frequency increases, especially during the last week of the simulation when the largest physical and numerical mixing occurs (Figure 7d).
There are a number of factors that could cause the difference between on-and offline estimates of numerical mixing, and many of them are difficult to quantify. First, we cannot say with certainty that  will yield identical results whether s or s′ is computed locally despite our analysis with a 1D upwind scheme. Another source of error is our use of lower-order accurate discretizations for offline analysis of the tracer variance budgets (see Appendix A). For example, ROMS has implemented several higher order internal time-stepping schemes for tracers since its inception (e.g., third-order Adams-Bashforth), whereas we use a first order centered finite difference for time derivatives in the s 2 and ′ 2 budgets. Likewise, we used volume-conserving fluxes in the calculation of the s 2 and ′ 2 boundary advection instead of reconstructing the MPDATA scheme offline. These discretization errors may be further compounded by our use of average files compared to snapshots.
The reason for this is the offline method is intended to focus on practicality. Recreating a model's complex numerical schemes offline may be more cumbersome than coding the online method into the source code and rerunning it if computational resources are not a concern. Other sources of error are the numerical diffusion of s 2 and ′ 2 through the lateral boundaries of the control volume or the generation of spurious convection. Quantification of numerical diffusion would require estimating a numerical diffusion coefficient and as we have shown, we cannot guarantee an offline calculation will be robust. The situation is similar regarding the quantification of spurious convection.
The online method should be used as the reference due to the lack of convergence and influence of discretization errors associated with the offline method. We can say that increasing the model output frequency will marginally improve the offline estimates, and that the ′ 2 budget will generally yield more accurate estimates of the numerical mixing than the s 2 budget, especially at low output frequency. In our case,  ′ 2 almost always overestimates  and qualitatively captures the temporal variability well. However, there is no guarantee this will be the case for other ocean models. Therefore, we cannot recommend the offline method for generic analysis of numerical mixing because it does not converge even at impractically high output frequencies and we are unable to identify the exact reasons for the lack of convergence. These issues may 10.1029/2022MS003380 23 of 30 be compounded for larger-scale ocean models. The proper output variables would have to be chosen before running the simulation and the output frequency will be decreased for larger-scale models. Future research could potentially improve the offline method by using higher-order numerical schemes for the tendency and advection terms, but this would have to be tested on a per-model basis and referenced against an online benchmark for validation.

Relating Numerical Mixing to |∇ H s| and Other Processes
As shown in Figure 5, numerical mixing may be over an order of magnitude larger than the physical mixing in regions with strong density gradients, generally corresponding to areas of larger |∇ H s|. Smolarkiewicz (1983) showed that after discretizing a one-dimensional advection equation with a first-order upwind scheme, the numerical mixing  is  where |u| is the magnitude of the constant horizontal velocity, Δx is the horizontal grid resolution, = Δ Δ is the Courant number with online time step Δt. After rearranging, it can be shown that the right-hand-side of Equation 26 is formulated such that the numerical mixing is equal to an implicit diffusion coefficient |u|Δx − u 2 Δt multiplied by the square of the horizontal salinity gradients. As Wang et al. (2021) notes, when the Courant number is less than one, the numerical mixing is approximately proportional to the square of the horizontal salinity gradients. Wang et al. (2021) applied this equation to a realistic simulation of the Changjiang River plume using two tracer advection schemes (MPDATA and Third-order Upstream-biased Horizontal Scheme) and suggested that this relationship holds true qualitatively (their Figures 8 and 9).
We further investigate the relationship between  and (∇ ) 2 = ( ) 2 + ( ) 2 in Figure 11, which shows weighted histograms of the absolute value of  and (∇ ) 2 in log 10 space, weighted by grid cell volume dV, for the coarse and fine simulations. We performed a weighted linear regression analysis to test the robustness of the two-dimensional form of Equation 26. There is a clear log-log relationship between  and (∇ ) 2 (r 2 = 0.55 for coarse simulation, r 2 = 0.39 for the fine simulation). To test for a power law dependence, we conducted an empirical linear fit of  and (∇ ) 2 in log 10 space, which slightly improved the fit (r 2 = 0.60 for coarse simulation, r 2 = 0.46 for the fine simulation). The relatively high r 2 suggest that the horizontal salinity gradients could be used to roughly approximate the numerical mixing, even for higher order and more complex advection schemes.
As (∇ ) 2 increases, more of the domain volume is concentrated at larger values of  . The relationship begins to taper off at (∇ ) 2 ∼  ( 10 −6 ) g 2 kg −2 m −2 , corresponding to salinity gradients associated with surface fronts and the pycnocline in the M/A river plume. These fronts have strong salinity gradients and numerical mixing, but they occupy a small portion of the water column, hence the decrease in grid cell volume. Weaker salinity gradients are associated with longer length and time scales and are correlated with weaker numerical mixing.
The differences between the coarse and fine simulations demonstrate that the fine simulation has weaker numerical mixing distributed at stronger salinity gradients than the coarse simulation. In other words, for a given salinity gradient, the fine simulation will experience less numerical mixing on average compared to the coarse simulation. The dividing line separating the positive and negative changes to dV between the simulations has a slope of nearly one, suggesting that the change in numerical mixing is proportional to (∇ ) 2 . The coefficient of determination is lower in the fine simulation because the numerical mixing depends not only on |∇ H s|, but the components of grid resolution dV, water velocities, and the model timestep. In the fine simulation, Δx and Δy decreased by a factor of five, but the time-averaged volume-integrated numerical mixing decreases by only 35%, suggesting a nonlinear relationship between numerical mixing, the square of the horizontal salinity gradients, and horizontal grid resolution. Although identifying the exact dynamical feedbacks that cause this non-linearity to arise is beyond the scope of this paper, an analysis (not shown) of the instability angle derived in Thomas et al. (2013) (their Equation 8) suggests that the fine simulation is more susceptible to symmetric and inertial instabilities than the coarse simulation. Coupled with the sharp changes of the velocity gradient tensor but insignificant changes to |∇ H s| in Figures 3 and 4, this is consistent with the idea that different dynamical processes emerge as the resolution is 10.1029/2022MS003380 24 of 30 increased, which complicates the relationship between horizontal resolution and numerical mixing. However, had |∇ H s| increased significantly, this would have increased the numerical mixing in the fine simulation.

Summary and Conclusions
We have studied physical and numerical mixing in a numerical model of the Texas-Louisiana (TXLA) continental shelf using a combination of on-and offline methods based on salinity variance. Physical mixing is defined as the dissipation of salinity variance associated with numerical closure schemes and numerical mixing Figure 11. Histograms of the horizontal salinity gradient squared (∇ ) 2 and absolute value of online numerical mixing | | weighted by grid cell volume dV for the coarse (a) and fine (b) simulations, and their differences (c). The thick dashed lines in (a) and (b) display weighted linear regression results for the approximate two-dimensional form of Equation 26 (black) and an empirical fit in log-log space (gray). The regressions and weighted coefficients of determination were obtained by subsampling the coarse grid model every three ξ, η points and the fine simulation every 15 points for the entire simulation period. The thick dashed gray lines in (c) indicate a slope of 1.

10.1029/2022MS003380
25 of 30 is defined as the mixing generated by the discretization of salinity advection. Salinity variance can be defined in terms of salinity squared s 2 and volume-mean salinity variance ′ 2 . Previous research (Burchard et al., 2019;MacCready et al., 2018; has shown the residuals of the s 2 and ′ 2 budgets can be used to estimate numerical mixing. However, the robustness of the offline method in realistic simulations has only been assessed qualitatively by Wang et al. (2021). The online method, which locally calculates numerical mixing  as the difference between the advected salinity squared and the square of the advected salinity divided by the model timestep (Burchard & Rennau, 2008).
 is used to evaluate the accuracy of the offline method.
We find the residuals of the s 2 and ′ 2 budgets,  2 and  ′ 2 , do not converge to the online method. This is true even as the model output frequency is increased to 10 min, an impractical output frequency for long-term realistic coastal models at this time. The s 2 budget at hourly output may miscalculate numerical mixing by over an order of magnitude. During the study period, 2 ≫ ′ 2 , which causes the resulting tendency and advection terms in the s 2 budget to be over an order of magnitude larger than the ′ 2 budget. We derive the s 2 budget in terms of volume-averaged salinity and salinity perturbation s′ to relate  2 and  ′ 2 and investigate the consequences of this scaling. We find the s 2 budget experiences larger truncation error compared to the ′ 2 budget when using identical numerical schemes.  2 begins to converge to  ′ 2 as output frequency increases but is still noisy at 10 min output. The scaling of s 2 and ′ 2 over the shelf is quite different than previous estuarine models (X. Li et al., 2018;L. Li et al., 2021;Warner et al., 2020), where ′ 2 may be similar in scale to s 2 since estuary domains include a transition from fresh river water to background coastal salinities.
The time-averaged  ′ 2 is 60% larger than the time-averaged  . Although  ′ 2 is much less sensitive to model output frequency, it does not converge to  for all times and cannot be considered robust. There are many sources of uncertainty that could contaminate the on-and offline methods. For example, a mismatch between our model's in-and external schemes for calculating time derivatives or additional sources of spurious mixing such as grid-scale noise in velocity.
Consequently, we cannot recommend the offline method for generic quantification of numerical mixing. Although Wang et al. (2021) qualitatively showed success in estimating numerical mixing offline, we find offline mixing estimates to be inaccurate, and only useful in providing a rough sense and scale of the numerical mixing. Despite this, the offline method qualitatively captures the larger signals in the temporal variability well and might be informative to researchers who are unable to use an online method. We also feel the approach might be useful in comparing different scenarios where the primary questions are about the relative magnitude of numerical mixing. This is particularly relevant for estuarine and coastal models already using a high model output frequency and spatial resolution, but less so for larger-scale models that employ a lower output frequency.
It is also clear that sources of uncertainty need to be evaluated on a per-model basis, so offline calculation cannot be generally recommended as a primary approach. Given the strong advection and weak physical mixing during summer in the nGoM, we think that our domain is particularly challenging for accurate offline calculation of numerical mixing, and is a good demonstration of the weaknesses of the offline approach. Other regions with relatively stronger physical mixing, for example, partially mixed estuaries, may be more suitable to the offline approach. However, a comprehensive understanding of when an offline approach may be feasible is outside the scope of the present study, and so for now we feel the offline approach must be treated as suspect for a particular region until demonstrated otherwise.
Regarding the online analysis, the numerical mixing remains significant relative to the physical mixing even for a submesoscale-resolving coastal ocean model. We find the volume-integrated numerical mixing comprises 57% of the bulk physical mixing. Instantaneously, the volume-integrated numerical mixing may exceed the physical mixing by almost half an order of magnitude, motivating us to use a two-way nested model with five times the native horizontal grid resolution to examine the sensitivity of numerical mixing to horizontal resolution. We find that the time-averaged volume-integrated numerical mixing decreases by approximately 35% in the fine simulation, suggesting a nonlinear relationship between horizontal resolution and numerical mixing. Building on the work of Wang et al. (2021), we use weighted property histograms to show that numerical mixing is approximately proportional to the square of the horizontal salinity gradients (∇ ) 2 . As horizontal resolution increases, this relationship weakens because newly resolved dynamical processes emerge, which are evident in histograms of the velocity gradient tensor and |∇ H s|. The surface s 2 flux surf n was calculated as where and are the evaporation and precipitation per unit area ("Evaporation" and "Rain" in ROMS output) computed online as part of the model forcing, and ρ w is the density of freshwater set to a constant value of 1,000 kg m −3 .
The volume-integrated horizontal diffusion Hdiff n was calculated as where κ H,i,j is the horizontal turbulent diffusivity computed offline defined in Equation 22. κ H,i,j was first linearly interpolated to the centers of cell edges to align the model grid after computing horizontal derivatives ∂ x and ∂ y of s 2 , which are calculated using a Jacobian to account for the change of the sea surface height with respect to time (Shchepetkin & McWilliams, 2005). Further details of the Jacobian calculation may be found in Thyng et al. (2022).
The resolved vertical mixing ( ) was calculated online so discretization errors do not bias comparisons with  . ( ) was calculated as where Vert. Mix is the resolved vertical mixing, denoted by "AKr" in ROMS syntax (not to be confused with AKs, the resolved vertical salinity diffusion coefficient). AKr was linearly interpolated to the cell centers of the vertical grid prior to integration. We note that offline calculations of ( ) , although not used in this study, tended to overestimate the online calculations.
The resolved horizontal mixing ( ) was computed offline because it is not currently programmed into the ROMS source code. It is also much smaller than other terms for our model, so offline discretization errors do not significantly impact  2 .
( ) was calculated as Horizontal gradients were calculated using a Jacobian but then interpolated to cell centers for both horizontal and vertical grids. Once calculated, the above terms may be substituted into Equation 8 to obtain  2 . The process for discretizing the volume-mean salinity variance budget is similar, but requires calculation of the volume-averaged salinity . was calculated as The volume-mean salinity variance ( ′ 2 ) may be calculated as where applicable into Equations A1, A2 and A4-A6. Recalculation of the physical mixing terms is not required because has no spatial gradients.

Data Availability Statement
All TXLA model output used in this study is publicly available at https://hafen.geos.tamu.edu/thredds/catalog/ catalog.html under the "TXLA ROMS nested model for SUNRISE/2010" subdirectory. The corresponding analysis code is available at https://doi.org/10.5281/zenodo.7566722. The calculations were performed in Python ver. 3.9 using the xarray (Hoyer et al., 2021), xgcm (Abernathey et al., 2022), and xroms (Thyng et al., 2022) packages. More information about specific packages used for analysis can be found in the code repository.