Application of the Forward Sensitivity Method to a GWCE-Based Shallow Water Model

The Forward Sensitivity Method (FSM) is applied to a GWCE-based shallow water model to analyze the sensitivity to the numerical parameter, G, that determines the balance between the wave and primitive forms of the continuity equation. Results show that the sensitivity to G calculated in the sensitivity evolution portion of the FSM is consistent with the actual sensitivity to G computed from multiple simulations using finite differences. The data assimilation step in the FSM is shown to be effective in selecting G that minimizes an objective function, in this case model errors based on sensitivities. Additionally, the FSM sensitivity results show 2∆x oscillations in the elevation and velocity fields develop when G is increased too high, suggesting the FSM may be an effective tool for determining the upper limit of G for real-world applications.


Introduction
The Generalized Wave Continuity Equation (GWCE) [1] is an extension of the Wave Continuity Equation (WCE) [2], which was introduced to eliminate the spurious oscillations that plagued finite element solutions of the primitive Shallow Water Equations (SWE).The GWCE contains a numerical parameter, G, that determines whether the GWCE tends towards a wave equation form (low G) or the primitive continuity equation (high G).It has undergone rigorous analytical studies, which have shown that the GWCE is consistent with the primitive continuity as long as the initial conditions satisfy continuity [1].However, Kinnmark went on to show that even if that condition is not satisfied exactly, then the solution remains robust as long as the numerical parameter, G, satisfies some minimal conditions, e.g., G > 0. Many other studies have shown the superior wave propagation characteristics of the GWCE, including non-aliased solutions for short waves (e.g., [3]) and low-dissipation for physical waves (e.g., [4]).In that G is a numerical parameter, akin to a penalty parameter commonly found in classic finite element methods, there have been numerous studies that sought to identify an "optimum" value of G (e.g., [5][6][7]).Additionally, it should go without saying that all numerical algorithms introduce conceptual errors (e.g., missing physics) and truncation errors into the solution; a goal of modeling is to minimize the adverse impact of those errors.A big contribution of the current manuscript is that it goes a step further than previous analyses because the methodology can be applied to nonlinear problems and because it opens the door for data assimilation, which is a widely-accepted practice of "tuning" a model to account for missing physics (e.g., subgrid scale processes).However, in the end, real-world applications over the last 30 years provide the truest test of the GWCE.For example, the resulting algorithm is employed in the production version of the widely-used ADCIRC code ( [8][9][10]), which has a long history of providing accurate, robust results in a wide variety of applications, including tide-and wind-driven circulation, hurricane storm surge and inundation, baroclinic transport, sediment transport and coastal dredging feasibility, and larval and oil spill conveyance settings ( [11][12][13][14][15][16][17]).
The WCE was first introduced by Lynch and Gray in 1979 [2]; in 1986, Kinnmark generalized the WCE to the GWCE by introducing a weighting factor, G, that is distinct from the bottom friction parameter, τ [1].Kolar et al. [5] found that G has a large effect on model results and that a value G > τ is necessary to minimize errors.Atkinson et al. [7] analyzed the wave propagation characteristics of the GWCE-based SWE, and they found that the GWCE-based system is nearly identical to the primitive SWE, with a quasi-bubble velocity approximation [18], for a specific G parameterization.The dispersion analysis results of [7] have guided the recent selection of G (cf. [11,16]), where spatially-variable parameter selection has been employed for diverse applications.However, specification of a value (or parameterization) for G is an on-going issue.
In general, techniques applied to analyze the GWCE-based system have been limited to linear analysis (or analysis of the linearized equations), e.g., dispersion and Fourier analysis.These classic techniques are also limited to constant bathymetry domains and interior nodes.Herein, the Forward Sensitivity Method (FSM) [19] is applied to analyze the 1D, GWCE-based SWE.In this analysis, both constant and non-constant bathymetry cases are analyzed.As mentioned in [19], the FSM builds on sensitivity function analysis (e.g., [20]) and includes an optimization component that allows observations to be used to correct the model.The FSM is a deterministic data assimilation strategy for correcting forecast errors when a deterministic model is used in the analysis.Forecast error is defined by the difference between the model solution and the given (noisy) observation that the model is supposed to capture in the first place.A model can be either perfect or imperfect.Recall that a solution of a dynamic model depends on: (a) initial conditions; (b) the values of parameters; and (c) the boundary conditions.Since these three factors control the evolution of the model solution, these are collectively called "control".The goal of FSM is to find corrections to the control so as to drive the forecast errors as close to zero as possible in the least squares sense.FSM was first reported in Lakshmivarahan and Lewis [19] and is closely related to the now classic adjoint sensitivity-based 4D VAR method [21].The FSM-based approach is quite general and can handle both linear and nonlinear models and can be used to correct the forecast errors due to all three components: initial conditions, boundary conditions and parameters.A comprehensive account of FSM and varied applications is given in Lakshmivarahan et al. [22].The method is applied to analyze a differential equation describing the air/sea interaction in [19].In contrast to dispersion analysis, which is limited in applicability (e.g., linear equations, interior nodes, constant bathymetry), the FSM can be applied to analyze the non-linear equations at all nodal locations within the domain.The FSM has the added capability of accounting for boundary conditions, whereas other methods look only at interior points.
While FSM is applicable to non-linear systems, the analysis in this manuscript is limited to the linear system, with the intent being to present the exploration of a new analysis tool for shallow water equation models.Application to the non-linear GWCE-based shallow water equations has been performed [23], and the results will be presented in a subsequent paper.As presented first, derivation of the equations for the evolution of the sensitivity functions follows [19].Then, the FSM sensitivities are analyzed for two domains, which is followed by a validation of the FSM sensitivities with a numerical analog sensitivity approach.Section 3 begins with the presentation of the methodology, based on [19], for computing parameter corrections, and concludes with applications for the linear sloping domain.Section 4 contains the results of a proof-of-concept sequential optimization.Subsequently, a comparison between FSM and dispersion analysis is presented in Section 5. Finally, conclusions are made based on the analysis herein.

Derivation of Sensitivity Equations
The 1D linear inviscid GWCE and momentum equation are given by Equations ( 1) and (2), respectively, where ζ is the water surface elevation, G is the weighting parameter in the GWCE, τ is the bottom friction term, h is bathymetry, u is depth-averaged velocity and g is the gravitational acceleration.Additionally, the subscripts denote partial derivatives, i.e., ζ tt is the second partial derivative of ζ with respect to time.Application of the continuous Galerkin finite element method, using constant grid spacing, and a finite difference time discretization results in Equations ( 3) and ( 4) for the GWCE and momentum equation, respectively, where, on an element basis, The system can be written symbolically, as: where the coefficient matrices A(G), B(G), C(G) ∈ R 2n×2n are square matrices with dimensions of twice the number of nodes, n, for each of the three time levels; the vectors of variables are c k+1 , c k , c k−1 ∈ R 2n ; and the forcing vector is f k+1 bc ∈ R n .The FSM allows calculation of the sensitivity to different aspects of the control, which includes initial and boundary conditions, as well as physical, empirical and numerical parameters.Herein, the focus is on the numerical parameter G.The sensitivity to G is the rate of change of the solution due to a change in G. Given the system described in Equation ( 5), the sensitivity is found by taking the derivative with respect to G, as shown in Equation (6).
Application of the product rule, the definition of the sensitivity of the solution to G at a given time as w k = ∂c k /∂G, and rearrangement yields: Note that the forcing vector is considered to be independent of G.According to Equation (7), the unknown sensitivity vector can be computed from the previous sensitivities and elevation and velocity fields, although c k+1 must be calculated before w k+1 .The three time-level scheme requires sets of sensitivity values at times k and k − 1. Results herein have cold start initial conditions, where the initial elevation and velocity fields are zero throughout the domain.As such, the initial conditions do not depend on G, and the initial conditions for the sensitivity to G are w −1 = w 0 = [0, . . ., 0] T ∈ R 2n .

Sensitivity Results for Tidal Problem on the Linear Sloping Domain
The parameters for the test case on the linear sloping domain are given in Table 1.The number of nodes (and, thus, the grid spacing) was varied, with 11, 21, 41 and 81 nodes, constant ∆x grids being employed.Additionally, the simulation duration and time step were variable.Finally, both explicit (α 1 = 0, α 2 = α 3 = 1/2) and implicit (α 1 = α 2 = α 3 = 1/3) versions of the code were assessed.The differences between the results from the explicit and implicit models (both flow variables and sensitivities to G) were immaterial over the stable range of G values, although the implicit α specification allows stability at lower G values than the explicit version (for this test case, the implicit model was stable at G values two orders of magnitude smaller than the lowest stable value using the explicit version).Simulations with each grid, using a G value of 0.001 s −1 were performed for a period of 3.0 days, with output recorded every 5.0 min for the last day.For the three coarsest grids, ∆t = 1.0 s, while a time step of 0.5 s was used for the 81-node grid.Nodal elevation and elevation sensitivity to G results at select locations in the domain are shown in Figure 1.The node number listed on each panel corresponds to the node number for the location in the 21-node grid.The first panel, labeled "Node 1", shows the specified (i.e., Dirichlet) elevation boundary time series.On each panel, there are four solid lines.Each line shows the temporal evolution of the water surface elevation at the specified location for a particular domain, corresponding to the 11-, 21-, 41-and 81-node domains.The four solid lines are overlain on one another because the time series at the boundary are equivalent.Additionally, each panel has four dashed lines, corresponding to the same grids as for the solid lines.The dashed lines for Node 1 (along the line y = 0) show that the elevation sensitivity to G is zero, which is due to the elevation boundary condition being independent of G.
The second panel, labeled "Node 3," shows results 4.0 km into the domain.For this test case, the elevation time series is independent of ∆x, as is evident by the indistinguishable solid lines.However, the magnitude of the elevation sensitivity to G is dependent on grid resolution, with the magnitude decreasing substantially with increased resolution.The reduction in sensitivity to G with increased grid resolution (i.e., smaller ∆x) suggests that a solution with only a limited dependence on G can be obtained for this domain if sufficient resolution is utilized.Additionally, the timing of the sensitivity is consistent for the different grids, with co-located zero sensitivity values, which are approximately 90 degrees out-of-phase from the zero elevation values.
Similar general trends hold for elevation and elevation sensitivity to G time series in the middle ("Node 11") and on the right side ("Node 21") of the domain.Again, the elevation time series are indistinguishable.However, the magnitude of the sensitivity to G is highly-dependent on ∆x.Additionally, the magnitude of the sensitivity to G, as well as the amount the elevation and elevation sensitivity time series are out-of-phase depend on the location in the domain, with the magnitude and the phase difference increasing with distance from the ocean boundary.Elevation (m) and elevation sensitivity to G (ms) for different simulations on the linear sloping domain, with each simulation using a different resolution grid (11-node, black; 21-node, red; 41-node, green; and 81-node, blue).The solid lines depict the elevation results, while the dashed lines show the temporal evolution of the sensitivity of the elevation to G. The node number listed in the title for each panel is the node number in the 21-node grid associated with a given location.
Figure 2 shows the velocity equivalents to the elevation results shown in Figure 1.At the ocean boundary, the velocity results are slightly dependent on ∆x.Additionally, it is noteworthy that there is a relatively large sensitivity to G at this location.The elevation value is specified at the boundary, so changes to velocity resulting from changes to G result in changes to the amount of mass entering and exiting the domain at the open boundary throughout the simulation.For the other locations in the domain, the velocity time series overlay one another.The velocity sensitivity to G at the ocean boundary shows that, regardless of grid resolution, the velocity is highly dependent on the choice of G, which has significant implications on global mass balance, as noted in [5].Throughout most of the domain, the phase shift of the velocity sensitivity to G is independent of ∆x, as was the case with the elevation results in Figure 1; the location denoted by "Node 3" is the aberration, as there is a phase shift for the velocity sensitivity for different grid resolutions.Furthermore, as with the elevations, the magnitude of the sensitivity to G decreases with increasing grid resolution, and the sensitivity is lower in magnitude the closer to the land boundary where velocity is specified.Velocity (m/s) and velocity sensitivity to G (m) for different simulations on the linear sloping domain, with each simulation using a different resolution grid (11-node, black; 21-node, red; 41-node, green; and 81-node, blue).The solid lines depict the velocity results, while the dashed lines show the temporal evolution of the sensitivity of the velocity to G. The node number listed in the title for each panel is the node number in the 21-node grid associated with a given location.
In order to assess the impact of different G values on the sensitivity of the solution to G, the implicit version of the code was used, with a ∆t of 5.0 s, for simulations with G values of 0.00001, 0.0001, 0.001, 0.01 and 0.1 s −1 .Additionally, the simulations were 10.0 days in duration.The sensitivity values over the last two days of the simulation, for Nodes 2-7, are shown in Figure 3.The gaps in Figure 3 correspond to times when the sensitivity value is below the minimum ordinate value (which is just greater than zero) on the plot, although generally, these instances correspond to times when the sensitivity is negative for the current set of simulations.
Increasing G results in a decrease in the magnitude of the peak sensitivity.To a lesser extent, increasing G changes the timing of the sensitivity.Specifically, when G is increased from 0.00001 to 0.0001 s −1 , there is a small decrease in the magnitude of the peak sensitivity and a shift in the timing, so the peak sensitivity occurs earlier.For the even-numbered nodes, these trends continue for subsequent increases in G, although the decrease in sensitivity magnitude is more prevalent than the shift in timing of the peak.In contrast to the results for the even-numbered nodes, the sensitivity results for the three highest G values at the odd-numbered nodes are not coincident in time.Specifically, at Node 3, the results for G values of 0.01 and 0.1 s −1 show a phase shift compared to the G value of 0.001 s −1 .
For the same set of simulations (varying G value), the velocity sensitivity to G follows the same general trends as the elevation sensitivity.The GWCE was introduced for CG finite element modeling to control spurious 2∆x oscillations present in solutions of the shallow water equations using the primitive continuity equation.Increasing G shifts the GWCE towards the primitive continuity equation.The elevation sensitivity results show 2∆x oscillations in the sensitivity to G for values of the numerical parameter of 0.01 s −1 and larger for this application on the linear sloping domain, suggesting those values result in the GWCE becoming "too primitive" for this test case.The decrease in the magnitude of the sensitivity to G as G increases is consistent with the formulation of the equations.Introduction of non-zero G values results in the primitive continuity portion of the GWCE contributing.Eventually, as G values are increased, the primitive continuity term becomes dominant, and further increases in G will have only minimal impacts on the solution.

Sensitivity Results for Tidal Problem over a Seamount
The simulation parameters for this second test case, using a seamount domain, are listed in Table 2; the values are similar to those used for the case in Section 2.2.The base, 31-node, seamount domain is shown in Figure 4.The length of the simulations was 5.0 days, and the time step was 5.0 s.The elevation sensitivity results (not shown), indicate that the general trends present from the linear sloping domain test case also apply for the seamount domain.Specifically, the magnitude of the sensitivity decreases with increasing G, and the peak sensitivities occur earlier in time for higher G values.Furthermore, node-to-node oscillations in the sensitivities occur for the higher G values in the set, with the 2∆x oscillations readily apparent for the highest G value, 0.1 s −1 .
The elevation sensitivity results from the four simulations with different G values are summarized by the left panel of Figure 5, which shows the peak elevation sensitivity to G, over the last day of the simulation, for each node in the domain for simulations with different G values (0.0001, 0.001, 0.01, 0.1 s −1 ).The general trend is for the peak elevation sensitivity to increase with distance from the ocean boundary.The results with G = 0.01 s −1 show 2∆x oscillations in the magnitude of the peak sensitivity for a substantial portion of the domain, which is indicative of the GWCE becoming "too primitive", even though the sign of the sensitivity does not follow the traditional 2∆x oscillation pattern that occurs for the highest G values.
The right panel of Figure 5 shows the peak velocity sensitivity to G over the last day of the simulation for each node in the seamount domain for the four simulations with different G values.The sensitivity is zero at the land boundary (Node 31); the peak velocity sensitivity increases from a minimum at the land boundary to a maximum over the seamount (Nodes 16-21), then decreases oceanward of the seamount.The results with a G value of 0.01 s −1 show short wavelength oscillations in the peak velocity sensitivity for the majority of the domain.In contrast, the results for the highest G value, 0.1 s −1 , do not show prevalent oscillations in the peak velocity sensitivity landward of the start of the rise of the seamount (Node 11).However, a smooth set of peak velocity sensitivity points is not a sufficient condition to conclude that the G value is below the "too primitive" threshold.Time series analysis of the velocity sensitivities for the highest G value reveals the node-to-node switching of signs on the sensitivities, i.e., the positive peak sensitivities for the odd-numbered nodes correspond to the same times as the maximum negative sensitivities for the even-numbered nodes, and vice versa.

Comparison of FSM and Numerical Analog Sensitivities
The FSM sensitivity results presented previously predict the changes in the solution (elevations or velocities) that result from a change in the numerical parameter, G.To verify the procedure for computing the FSM sensitivity, the FSM sensitivity is compared to a numerical analog sensitivity.The numerical analog sensitivity is computed using finite differences.In particular, by comparing the results from two simulations with different G values, finite difference approximation of the sensitivity to G can be calculated using Equation (8), where the subscripts 1 and 2 refer to the two different simulations.
A comparison of the FSM and numerical analog elevation sensitivities to G for Node 11 in the linear sloping domain is shown in Figure 6.The left panel in Figure 6 is a comparison of the FSM sensitivity to G (black line), for a simulation with a G value of 0.001 s −1 , to the numerical analog sensitivity (red line) calculated using results from simulations with G values of 0.001 and 0.003 s −1 (∆G = 0.002).The evolutions of the sensitivities have the same shape, although the magnitude of numerical analog is significantly lower than the FSM sensitivity.The right panel of Figure 6 is a comparison of the FSM sensitivity to G to the forward numerical analog with a smaller difference in G values, ∆G = 0.0001 s −1 .It is readily apparent that decreasing the difference in G used to compute the numerical analog reduces the difference between the FSM and numerical analog sensitivities, although the numerical analog sensitivity is still smaller in magnitude than the FSM sensitivity.However, this underprediction by the numerical analog is directly related to the choice of G values used to compute the numerical analog.In this case, the second value of G used to generate the numerical analog is larger than the value of G for which the sensitivity is desired, which will be referred to as a forward numerical analog (because G 2 > G 1 , meaning ∆G is positive).The sensitivity to G decreases with increasing G (c.f., Figure 3), so the forward numerical analog to G is generally lower than the FSM sensitivity for a simulation with G = G 1 .Furthermore, for forward numerical analogs, increasing ∆G increases the underprediction.In contrast, use of a similar backward numerical analog would show that the numerical analog sensitivity is slightly greater in magnitude than the FSM sensitivity.
The results presented above show that the sensitivities computed using the FSM are consistent with the sensitivities calculated using the numerical analog as ∆G goes to zero, and the comparison confirms that the behavior predicted by the FSM actually occurs in the solution as G varies.As such, the FSM presents an opportunity to perform data assimilation, explored in Section 3, based on errors between observations and results from a simulation with a given value of G, although one could alternatively use a numerical analog approach to compute sensitivities for use in the data assimilation step.The equivalence of the FSM and numerical analog sensitivities gives rise to the following question: what is the benefit of FSM over a reasonably simple and straightforward numerical finite difference calculation?In the case of the constant G simulations presented herein, the two methods would require similar computational effort.However, for multi-parameter estimation, as is required for spatiallyand temporally-variable G specification, where sensitivities to p variables is necessary, p + 1 ADCIRC simulations would be necessary to compute the p numerical analog sensitivities (with one base run and p simulations with a small change in each parameter), whereas the FSM sensitivities to p parameters can be computed during an individual simulation.

Data Assimilation Approach
The second component of the FSM is the data assimilation step to correct G using the sensitivities and computed model errors.As presented in [19], using a first-order approach, where only the first term in the Taylor series expansion is retained, the error is equal to the product of the sensitivity and the correction, as given by Equation ( 9), e(x, t, G) = z(x, t) − c(x, t, G) ≈ (∆G)w(x, t, G) where for spatial location x, time t and numerical parameter value G, e(x, t, G) is the simulation error, z(x, t) is the observation value, c(x, t, G) is the model result, ∆G is the correction to the numerical parameter and w(x, t, G) is the sensitivity to G.
The correction can be computed in a variety of ways.The simplest correction uses an observation at one point in space, x j , at one time, t k , along with the model results for the same location in space and time.The correction, ∆G, to the value used for the simulation, G 0 , based on this one observation is shown in Equation (10).
Least-squares minimization is a more sophisticated approach that allows for the use of multiple observations in space or time.For the results herein, least-squares minimization will be applied on a nodal basis.In other words, the observations and model results for a given point in space, over a range of time, will be used to compute a least-squares correction to G.This is analogous to the real-world situation where a buoy collects a time series of water surface elevation data at a fixed location in the domain.Conversely, least-squares minimization could be applied on a temporal basis where errors throughout the domain, at a given time, are used to generate a correction to G.
The least-squares correction, based on results and observations for node j using nrecs values in time, requires the vector of sensitivities H j and the error vector e j .
The optimal least-squares correction, adapted from [19] for a scalar parameter, is given by Equation ( 13).
(∆G) j = H j T e j H j T H j (13) The optimal least-squares correction is a standard result that is presented in [21], which provides additional detail about the origins of the analysis technique.

Correction to CG Results on the Linear Sloping Domain
In this section, "observations" will be taken from model results generated using the 2D CG version of ADCIRC on a rectangular grid that is uniform in the y-direction.The 2D code was run implicitly with the same parameters as the 1D code, and the 2D domain consists of 11 nodes in the y-direction for each of the 21 nodes in the x-direction for the linear sloping domain.Results for the sixth line of nodes (the centerline) from a simulation with a constant G value of 0.001 s −1 are used as the observations.Furthermore, the x-component of the velocity from the 2D model is used as the velocity observation; the y-component of the velocity is ignored, but is generally several orders of magnitude less than the x-component (and close to zero).
The purpose of the data assimilation step in the FSM is to reduce model error.Therefore, before delving into the calculations of the corrections, it is informative to analyze model error for a range of G values.The error metric is the temporal mean of the root mean square error in space, denoted as RMSE x .The equation for the temporal mean of the root mean square elevation error in space is shown in Equation (14).
The elevation error results are shown in the left panel of Figure 7, while the velocity error results are shown in the right panel.For the G values used in the implicit 1D simulations, the minimum elevation and velocity errors are achieved when approximately the same G value is used in the 1D model as was used in the 2D model (to create the observations).The value of G that minimizes the error (i.e., G ≈ 0.001 s −1 ) is the value of G that should be revealed using the data assimilation step of the FSM.When available, elevation data are often in the form of time series at discrete location.Thus, herein, a time series of elevation observations will be used to calculate errors, and the correction will be computed using those errors and the corresponding time series of sensitivity values for a given node, as per Equation (13).
The correction varies for a given run depend on which node is used to calculate the correction.For example, the corrections based on the results from the simulation with G = 0.0001 s −1 are shown in Figure 8.The error and the sensitivity of the elevation results to G are both zero at the left boundary node, so the correction is not computed at that location (Node 1); rather, the correction is set to zero for plotting purposes.The correction can be calculated for each of the other nodes in the domain, and Figure 8 shows the correction to be just slightly greater than 0.0001 s −1 for each of the nodes.However, we know the optimal correction is close to 0.0009 s −1 , based on the values of G used for the runs to generate the model and observation results.The discrepancy between the computed least-squares correction and the optimal correction (which would result in the new value of G being the one that minimizes the model error) is a result of the variation in the sensitivity with G, as well as the fact that only the first order terms are kept in the Taylor series development of the correction equation.The sensitivity to G is much greater when G = 0.0001 s −1 than when G = 0.001 s −1 .Because the correction varies inversely with the sensitivity, the correction calculated using the sensitivity from the run with G = 0.0001 s −1 is, expectedly, low.
In order to show how the correction varies with G, the maximum, minimum and mean of the nodal corrections were calculated.Referring back to Figure 8, which shows a set of nodal least-squares corrections for the simulation with G = 0.0001 s −1 , the maximum correction is from Node 2, ∆G = 0.000168 s −1 , while the minimum comes from Node 9 (∆G = 0.000107 s −1 ).The mean correction is the arithmetic mean of the nodal corrections for Nodes 2-21.The results are shown in Figure 9.As seen in Figure 9, for a given simulation, the maximum, minimum and mean corrections are similar for simulations with G values less than 0.005 s −1 .Thus, regardless of location in the domain, the correction is similar, as was the case for the set of corrections shown in Figure 8.Interestingly, for the simulations with G values less than 0.001 s −1 , the corrections are larger for the simulations with G values closer to the target value, which seems counterintuitive.However, as mentioned previously, the large variation in sensitivity with G causes under-corrections for simulations with low G values.For the simulations with G values of 0.002 and 0.005 s −1 , the corrections are consistent and negative, as expected.However, the mean correction is larger in magnitude than the value of G used for the simulation.For example, the mean correction for the simulation with G = 0.005 s −1 is ∆G = −0.289s −1 .In contrast to the corrections from simulations with G values less than the optimal value, simulations with G values greater than the optimal value have corrections that are too large in magnitude.
Furthermore, for G values of 0.01 s −1 and above, some of the corrections are positive (indicated by the black dots on the top left panel in Figure 9), which is opposite in sign from the mean corrections.The presence of positive and negative corrections for the same simulation is a result of the GWCE becoming "too primitive".The initial appearance (lowest G value that experiences oscillations) corresponds to the G threshold above which spurious oscillations are generated.When the solution becomes "too primitive", the sensitivities start to become irregular.Rather than being similar from one node to the next, the sensitivities for successive nodes are opposite in sign or have varying magnitudes of the same sign.This transition from a normal pattern of sensitivities to an irregular one produces the aberrant correction results.
The difference in results for two model simulations is given by Equation (15).
For G values just greater than G = 0.001 s −1 , the sensitivities are similar from node to node.Therefore, the errors for model simulations with small deviations from the target value used in these studies will be similar between nodes, as long as G is not increased too much.When G = 0.01 s −1 , there are 2∆x oscillations in the sensitivities.Thus, at some G value between 0.001 and 0.01 s −1 , 2∆x oscillations begin to develop in the error values as a result of the oscillations in the sensitivities.
By computing the numerical analog using the target value as one of the G values for the simulation, the result is the average sensitivity over the span of G values.This average sensitivity can be compared to the FSM sensitivity, which gives the instantaneous sensitivity value.Figure 10 is a comparison of the numerical analog sensitivity between G values of 0.001 and 0.1 s −1 and FSM sensitivity for G = 0.1 s −1 for the 11th and 12th nodes in the linear sloping domain.It is readily apparent that the FSM sensitivity results are opposite in sign for the two nodes.However, the numerical analog sensitivity results are similar for the two nodes.The notable difference is the magnitude of the numerical analog sensitivities is larger for Node 12 than Node 11, which implies there is more error for results at Node 12 than Node 11.In this case, the ∆G value used to compute the numerical analog is −0.099 s −1 .Therefore, when the numerical analog sensitivity is positive, the error is negative, and vice versa.It should also be noted that the FSM sensitivities are, generally, in-phase with the numerical analog sensitivities for Node 12, whereas the two sets of sensitivities are out-of-phase for Node 11.Therefore, an additional increase in G away from G = 0.001 s −1 will cause increases in the magnitude of the error at Node 12 and decreases in the magnitude of the error at Node 11.
As mentioned previously, the occurrences of positive numerical analog sensitivities in Figure 10 (e.g., the peak values occurring approximately 8.0, 8.5, 9.0, 9.5 and 10.0 days into the simulation) correspond to times of negative model error (compared to the simulation with G = 0.001 s −1 ), and vice versa.Therefore, 9.0 days into the simulation, the error is negative at Nodes 11 and 12.For Node 11, the numerical analog and FSM sensitivities are out-of-phase, which means that, generally, when the FSM sensitivity is positive/negative, the error is positive/negative (numerical analog sensitivity is negative/positive).Subsequently, the correction to G will be positive, which is the wrong direction.In contrast, the numerical analog and FSM sensitivities are generally in-phase at Node 12, which results in the correction to G being negative, because the product of the error and sensitivity vectors is negative.The corrections to G for each of the nodes are shown in Figure 11.As expected, based on Figure 10, the correction produced using results for Node 11 is positive, while the correction generated using results for Node 12 is negative.

Sequential Optimization
In the previous section, corrections to G were calculated based on model errors and sensitivities to G.In this section, the correction, ∆G, is added to the previous G value to determine the next G value.This process is continued until the new correction is below a certain threshold, which signifies that the optimization process has converged at the target value.
The linear sloping domain is used for this proof-of-concept application, along with the explicit version of the code.The simulation parameters are the same as those used previously, with the exception that the run is only 5.0 days long, and corrections are generated using the results from the last day of the simulation.The observations are the elevation results along the centerline of the 2D ADCIRC simulation with G = 0.001 s −1 .The correction, ∆G, is the mean of the nodal corrections using the elevation results to compute the errors.
The initial G value for this exercise is 0.0005 s −1 , and the convergence threshold for ∆G was set at 1.0 × 10 −10 s −1 .As expected, specification of an initial value that is less than the target value resulted in each correction being in the appropriate direction (positive), with less than the optimal magnitude, as shown in Table 3.The target value is the value at which the sequential optimization finishes, 8.98 × 10 −4 s −1 (however, it is close to the value of G used in the 2D model to create the observations, but in practice, this would not be known).The ratio of the correction to the optimal correction is notable; as the G value approaches the target value, the correction approaches the optimal correction.This is logical because as the difference between the current and target G values goes to zero, the instantaneous sensitivity to G given by the FSM gets closer to the average sensitivity over the range.Table 3. Sequential optimization of G for the linear sloping domain compared to the 2D CG simulation with G = 0.001s −1 .The units for each of the columns, except for the fourth column, are s −1 .The ratios in the fourth column are dimensionless.Additionally, tests were performed with an initial G value greater than the target value.As expected, specification of a value larger than the target, but still below the primitive threshold, results in an over-correction in the first step.For the explicit code with an initial value of 1.20 × 10 −3 s −1 , the mean correction for the first step is ∆G = −3.94× 10 −4 s −1 , resulting in a new G value of 8.06 × 10 −4 .From there, the corrections bring the G value up to the target value.However, if the initial specification is significantly higher than the target, the over-correction can result in negative G value.For instance, the mean correction with G = 2.00 × 10 −3 s −1 is −2.38 × 10 −3 s −1 , which is larger than the previous G value.Thus, in practice, constraints on G would have to be put into the optimization algorithm.

Comparison of FSM to Dispersion Analysis
Kolar et al. [5] performed a dispersion analysis of the 1D shallow water equations using the GWCE for the Bight of Abaco, Bahamas.Kolar et al. noted that spurious 2∆x oscillations do not occur if the dispersion curve is monotonic.In their paper, they delineated the frequency for the M 6 tide [5] (p.536) and found that the monotonic dispersion relations for this frequency exist as long as G does not exceed 0.075 s −1 .The frequency of the M 2 tide is one third the frequency of the M 6 tide, so G must be less than approximately 0.3 s −1 to ensure the solution remains free of spurious, short-wavelength oscillations for the M 2 frequency.
The dispersion analysis performed in Kolar et al. [5] used a bathymetry value of 2.0 m, an element size of 2700 m and a bottom friction value of 0.01 s −1 .For this study, these parameters were also used in 1D simulations with a flat bottom domain consisting of 21 nodes.The time step for the 1D simulations was 5.0 s.Larger time steps result in differences in the calculated sensitivities, whereas the sensitivities were consistent between simulations with time steps of 2.5 and 5.0 s.It should also be noted that dispersion analysis is restricted to interior nodes.The 1D simulations herein using the FSM include boundaries that are treated as stated previously (specified elevation on the left, zero velocity on the right).
For the 1D simulations using the M 6 tide, the dispersion analysis predicts spurious oscillations for G values greater than 0.075 s −1 .The elevation FSM sensitivity results are free of 2∆x oscillations with G = 0.01 s −1 .With G = 0.03 s −1 , the sensitivity results show 2∆x oscillations for the first four elements in from the left boundary.However, the interior of the domain is not impacted.Further increase of G results in oscillations in a greater percentage of the domain.
Using the M 2 tide, the dispersion analysis predicts spurious oscillations for G values greater than or equal to approximately 0.3 s −1 .Again, oscillations in the FSM elevation sensitivities to G do not occur with G = 0.01 s −1 and occur only near the ocean boundary with G = 0.03 s −1 .Similar to the case for the M 6 forcing, the oscillations become more prevalent as G increases, although the M 2 forced simulations generally have less prominent 2∆x noise than the simulations with the M 6 forcing.This is consistent with the results suggested by the dispersion analysis.With a G value of 0.1 s −1 , the entire domain experiences 2∆x oscillations.
The FSM sensitivities and dispersion analysis do not produce exactly the same values of G for the onset of 2∆x oscillations in the solution.Given the underlying differences in the analysis techniques (e.g., dispersion analysis is confined to interior nodes and continuous time), this is not an entirely surprising result.However, the similarity between the results for the two techniques points to FSM being a useful tool in the analysis of problems where dispersion analysis is not valid (e.g., non-linear equations, non-constant bathymetry, etc.).

Conclusions
The FSM was successfully applied to the linearized, 1D version of ADCIRC with constant G.The FSM is useful in determining the sensitivity, both in space and time, to G.In particular, the sensitivity of the elevation and velocity fields to changes in G varies greatly with G.The sensitivity is much greater at low values of G than at higher values, where the GWCE effectively approaches the primitive continuity equation.Additionally, as G increases, the sensitivities from the FSM show the 2∆x oscillations that plague the continuous Galerkin finite element solution when the primitive continuity equation is used instead of the GWCE.Furthermore, the maximum G threshold, above which the GWCE becomes "too primitive" and results in the generation of spurious 2∆x oscillations, can be identified through analysis of the FSM sensitivities.In that sense, FSM can be used as a tool like dispersion analysis to predict the folding of dispersion relations, with the advantage of being applicable to complex, real-world problems.
The corrections, ∆G, calculated in the data assimilation step of the FSM are intrinsically tied to the sensitivities.The change in sensitivity over the range of possible G values makes direct estimation of the optimal correction difficult using first-order methods.At high G values, the corrections are also hindered by the 2∆x oscillations in the sensitivities.However, sequential optimization should be possible as long as care is taken in the specification of the starting point for optimization.Specifically, use of a low initial value is optimal because the corrections are more stable, compared to higher values of G.
While this analysis was limited to the linearized, 1D SWE, the FSM has potential use in more complex systems.Additionally, while the analysis is focused on the sensitivity of the system to G, the method can be adapted to analyze other parameters of the model.

Figure 1 .
Figure1.Elevation (m) and elevation sensitivity to G (ms) for different simulations on the linear sloping domain, with each simulation using a different resolution grid (11-node, black; 21-node, red; 41-node, green; and 81-node, blue).The solid lines depict the elevation results, while the dashed lines show the temporal evolution of the sensitivity of the elevation to G. The node number listed in the title for each panel is the node number in the 21-node grid associated with a given location.

Figure 2 .
Figure2.Velocity (m/s) and velocity sensitivity to G (m) for different simulations on the linear sloping domain, with each simulation using a different resolution grid (11-node, black; 21-node, red; 41-node, green; and 81-node, blue).The solid lines depict the velocity results, while the dashed lines show the temporal evolution of the sensitivity of the velocity to G. The node number listed in the title for each panel is the node number in the 21-node grid associated with a given location.

Figure 3 .
Figure 3. Elevation sensitivity to G, for simulations with different G values, at different locations on the 21-node linear sloping domain.The five different lines on each plot correspond to the five simulations, each with a different value of G as indicated in the legend below the figures.(a) Node 2; (b) Node 3; (c) Node 4; (d) Node 5; (e) Node 6; (f) Node 7.

Figure 4 .
Figure 4. Bathymetry and node locations for the seamount domain.

Figure 5 .
Figure 5. Peak elevation (left panel) and velocity (right panel) sensitivity to G for implicit runs on the seamount domain.Each dot denotes the peak elevation sensitivity value, for a given node, over the last day of a five day simulation.The color of the dot is based on the G value shown in the legend.

Figure 6 .
Figure 6.Comparison of forward sensitivity method and numerical analog elevation sensitivity results for the last two days of an explicit model simulation on the linear sloping domain.In each panel, the black line depicts the temporal evolution to the FSM elevation sensitivity to G for a simulation with G = 0.001 s −1 .The red line shows the time series of the numerical analog sensitivities.In the left panel, the two G values used for the simulations were 0.001 and 0.003 s −1 , while G values of 0.001 and 0.0011 s −1 were used to generate the numerical analog for the right panel.(a) ∆G = 0.002 s −1 ; (b) ∆G = 0.0001 s −1 .

Figure 7 .
Figure 7. RMSE x (ζ) (left panel) and RMSE x (u) (right panel) results, over the last two days of 10.0-day simulation, for simulations with different values of G on the linear sloping domain, using observations from the 2D CG ADCIRC model with G = 0.001 s −1 .

s 1 Figure 8 .
Figure 8. Nodal least-squares corrections, ∆G, for each node based on results for an implicit simulation with G = 0.0001 s −1 using output from the last two days of the 10.0-day simulation on the linear sloping domain.Observations are from the 2D CG ADCIRC code with G = 0.001 s −1 .

Figure 9 .
Figure 9. Maximum, minimum and mean nodal least-squares correction, ∆G, for simulations over a range of G values, using output from the last two days of the 10.0-day simulation on the linear sloping domain.The magnitude of the correction is shown as the ordinate, while the color of the dot corresponds to the sign of the correction: positive corrections are shown in black, and negative corrections are shown in red.Observations are from the 2D CG ADCIRC code with G = 0.001 s −1 .(a) Maximum; (b) minimum; (c) mean.

Table 1 .
Parameters for the linear sloping domain test case.

Table 2 .
Parameters for the seamount domain test case.