Definition and solution of a stochastic inverse problem for the Manning's n parameter field in hydrodynamic models.

The uncertainty in spatially heterogeneous Manning's n fields is quantified using a novel formulation and numerical solution of stochastic inverse problems for physics-based models. The uncertainty is quantified in terms of a probability measure and the physics-based model considered here is the state-of-the-art ADCIRC model although the presented methodology applies to other hydrodynamic models. An accessible overview of the formulation and solution of the stochastic inverse problem in a mathematically rigorous framework based on measure theory is presented. Technical details that arise in practice by applying the framework to determine the Manning's n parameter field in a shallow water equation model used for coastal hydrodynamics are presented and an efficient computational algorithm and open source software package are developed. A new notion of "condition" for the stochastic inverse problem is defined and analyzed as it relates to the computation of probabilities. This notion of condition is investigated to determine effective output quantities of interest of maximum water elevations to use for the inverse problem for the Manning's n parameter and the effect on model predictions is analyzed.

1. Introduction

In shallow water hydrodynamics, especially in coastal regions, channels and overland flows, bottom stress induced by surface roughness heavily influences flow behavior. Bottom stress is typically parameterized using the Gauckler-Manning-Strickler formula which relates average cross-sectional velocity, hydraulic radius, and the channel bed slope [1, 2, 3, 4, 5] and includes a coefficient that depends on the specific characteristics of the bottom surface according to standard definitions. In this paper, we use Manning’s formula to describe bottom stress and focus on the determination of the so-called Manning’s n coefficient, which is a spatially dependent field variable. Additionally, we investigate the indirect determination of the Manning’s n coefficient from observations on the state of water elevation and momentum, where Manning’s n coefficient is connected to wave motion through a deterministic shallow water model.

This is an example of an inverse problem whose solution provides the basis for scientific inferences for a deterministic, physics-based model. The goal is to predict unobserved behavior of a system using the mathematical model fed by input data and parameters that characterize the physical properties of a given state of the system. However, it is often infeasible or prohibitively expensive to experimentally observe all of the important input quantities that characterize the physical properties of a given state. Rather, the available data tends to be on observable aspects of the state of the system itself. Moreover, the set of observable quantities is generally different than the set of quantities to be predicted because of physical limitations in what, where, and when data can be collected.

In Fig. 1, we illustrate the abstract process of solving an inverse problem for a physics-based model in order to provide the input for scientific inferences based on model predictions. The physics-based model is given as a set of differential equations whose solution produces a map between the space of input data and parameters to the solution space. A set of functionals form a map from the solution space to the space of observable quantities, while another set of functionals map the solution space to the space of quantities we wish to predict. These are illustrated by the red arrows. Note that the observation and prediction functionals can be physically different transformations or they can be the same set of transformations but evaluated at different times and locations. To solve the prediction problem, we first solve the inverse problem for the map from the space of input data and parameters to the observation space obtained by composing the observable functionals with the solution of the physics-based model. We then solve the model “forward” with those data and parameter values and apply the prediction functionals to the solution(s) to produce the prediction. This flow is illustrated with black arrows.

Figure 1

The input space of data and parameters are mapped by the physics-based model to the solution space. Generally, different sets of functionals map solutions to quantities we can observe and to the quantities we wish to predict. Parameters and data for the model are obtained from observable data by solving an inverse problem for the model+observable functionals. Information on parameter space is then used to make predictions by “forward” computation.

The determination of Manning’s n for a model of coastal hydrodynamics provides a concrete example. Manning’s n characterizes the momentum loss due to bottom friction that affects bottom stress in a hydrodynamics model in Manning’s formula for drag. The Gauckler-Manning-Strickler formula was originally developed for open-channel flow in a fully-turbulent fairly regular unvegetated domain [1, 2, 3, 4, 5, 6]. For open-channel flow that meets these assumptions Manning’s coefficient may be determined indirectly from field or laboratory measurements of flow characteristics or from a sampling of the diameters of roughness elements [7, 8, 9, 10]. Arcement and Schneider [11], along with Chow [12], and Barnes [8] have provided guides for selecting Manning’s coefficient that extend the original domain of application to flood plains, vegetated water ways, and channels with irregular boundaries. These guides supply tables [12], photographic guides [8], and a combination thereof [11]. Furthermore, Manning’s formula, although extremely useful, is not the correct model of momentum loss due to vegetation in natural waterways as Manning’s formula models bottom friction not the complex form drag of flexible vegetation [13, 14, 15]. However, complex hydrological models require some representation of momentum loss due to a combination of bottom friction, vegetation, bed-forms, and the porous media-like structures that occur in coastal estuaries [16, 17, 18]. Detailed field measurements can be used to determine Manning’s coefficient through hydraulic calibration for a specific geographical location [19].

Unfortunately, it would be extremely cost prohibitive to use field measurements to estimate Manning’s n for these conditions at a fine detail over a large physical domain. One set of observation data that is available in coastal hydrodynamics is the maximum water surge heights at various fixed observation stations. We aspire to leverage the land cover and classification data collected by the NCLD and CCAP projects to estimate Manning’s coefficient for these types of computational modeling problems [20, 21, 14, 22]. This sets up the problem of determining the Manning’s n values from the maximum surge height data by solving the inverse problem for the shallow water hydrodynamics model. A related prediction problem is then to use the computed Manning’s n to predict inundation at certain critical points in the domain, e.g. near physical barriers intended to reduce the effects of flooding.

The inverse problem for prediction is complicated by two issues. First, the map from the data and parameter space to the observable space generally reduces the dimension, i.e. it is a “many-to-few” map. This means that the inverse problem has solutions, i.e. there is a set of possible data and parameters values corresponding to each observation data point. Second, all of the available data is subject to natural variation as well as experimental/observational error so the solutions of the inverse problem for parameter determination and the forward prediction problem are described in terms of probability measures. Both of these issues can be addressed directly by measure theory, which provides a very natural framework for the formulation, solution, and numerical approximation of the inverse problem for scientific inference [23, 24, 25, 26].

Measure theory is fundamentally based on the notion of an inverse of a map, and so the measure theoretic description of the stochastic inverse problem given below is universal in nature. There are other ways to formulate inverse problems for determining parameters in models. For example, the most common deterministic approach alters the inverse problem through “regularization” to obtain a new problem that has a unique solution. Numerically, the solution is often formulated as the solution of an optimization problem that is then solved approximately [27]. There are other techniques for solving stochastic inverse problems, e.g. Bayesian approaches [28, 29, 30] and Kalman filtering [31]. However, outside of our solution approach, we do not know of any stochastic solution methods that clearly adhere to the structure of the inverse problem, i.e. set-valued solutions, of the original model. It may be that they are in fact either solving a different inverse problem (as is the case when a parameterized distribution is assumed for the solution) or have altered the original problem (Bayesian approaches often implicitly involve regularization).

In the following, we let Q denote the map formed by the composition of the observable functions with the solution operator of the model, Λ denote the domain of input values for Q, and 𝒟 ≔ Q(Λ) denote the range. This Q should not be confused with the quantities Q and Q in the shallow water equations.

2. Description of the model and the Manning’s n parameter field

Assuming hydrostatic pressure and domains with large horizontal length scales relative to vertical length scales, the shallow water equations (SWEs) are used to model the hydrodynamic system. The SWEs can be derived by depth integration of the incompressible Navier-Stokes equations or directly from the integral form of the fundamental equations of fluid dynamics, resulting in a coupled system of equations consisting of a first-order hyperbolic continuity equation for water elevation and momentum equations for horizontal depth-averaged velocities.

Let ζ denote the free surface elevation relative to the geoid, then the continuity equation is given by and the momentum equations for depth-averaged x- and y-velocities (denoted by U and U, respectively) are given by Here, h is bathymetric depth relative to the geoid, H = ζ + h is the total water column height, f is the Coriolis parameter, P is the atmospheric pressure at the free surface, ρ0 is the reference density of water, and η is the Newtonian equilibrium tide potential. Letting x = x or x = y, then Q = U is the flux per unit width in the x direction, M is the vertically-integrated lateral stress gradient, D is the momentum dispersion, and B is the vertically-integrated baroclinic pressure gradient. Finally, τ are the imposed surface stresses and τ are the bottom stress components.

In regions where water depths are quite shallow; e.g., low-lying coastal areas, tidal flats, channels, etc., the bottom stress components τ are more sensitive to Manning’s coefficent (Eq. 2.2). These are defined using a linear or quadratic drag law through the coefficient K, We utilize a quadratic drag law, so K = C |U|. The constant C is determined using the Manning’s n formulation, where n denotes the Manning’s n coefficient of roughness. Note that for H sufficiently small C is set to a fixed constant to prevent division by zero.

3. Numerical solution of the hydrodynamic model and quantities of interest

We use ADCIRC to solve a hydrodynamic model that depends on 16 parameter fields, including bottom stress, bathymetric depth, and surface stress, using finite element methods on unstructured triangular meshes discretizing the spatial domains and using finite difference schemes in time. In ADCIRC, the continuity equation, Eq. (2.1), is replaced by a second-order, hyperbolic generalized wave continuity equation (GWCE) [32, 33] to reduce spurious oscillations that can occur in numerical solution of the form of the SWEs discussed above [34]. The steps for obtaining the GWCE are: (1) a multiple, τ0 ≥ 0, of the continuity equation is added to the time derivative of the continuity equation, (2) bathymetric depth is assumed constant in time, i.e., ∂H/∂t = ∂ζ/∂t, and (3) substitute the momentum equations into the continuity equation. The interested reader can find all details of these equations, their derivation, and numerical implementation in ADCIRC in [35].

The ADCIRC model has undergone extensive verification and validation, e.g., using hindcast studies with data from hurricanes over the last five decades including including Hurricanes Betsy (1965), Ivan (2004), Dennis (2004), Katrina (2005), Rita (2005), Gustav (2008), and Ike (2008) [36, 37, 16, 17, 38]. ADCIRC has also been used to generate tidal constituent databases for the eastern coast of the conterminous U.S. [39].

The verification and validation studies of ADCIRC have consistently shown that the maximum water heights obtained from the model are in excellent agreement with those observed in the field by fixed observation stations. Therefore, in this work, we consider the quantities of interest as the maximum heights at a given set of locations.

4. Formulation and solution of the stochastic inverse problem

4.1. Measure theory and inverse problems

Probability is introduced when we assume stochastic models for the natural variations and experimental error in data. However, this also introduces measure theory because rigorous probability is described in terms of measure theory [40]. Measure theory is based on the concept of inverse problems, hence a rigorous probabilistic description of the inverse problems for parameter determination in the presence of stochastic variation has a very specific mathematical formulation. There are features of this formulation that must inform any approximate solution technique for the original stochastic inverse problem for a deterministic model. There are other ways to formulate inverse problems for parameters in a model and solution techniques that involve altering the original stochastic inverse problem. We give brief descriptions of examples below at appropriate places. To our knowledge, a systematic comparison of different inverse problems for parameter determination has not been carried out.

Below, we summarize the ingredients of measure theory briefly, and then express the relevant results about inverse problems and probability measures in terms of iterated integrals in several dimensions. The expression of probability measures in terms of integrals is key to all stochastic or deterministic approximate solution methods. In this way, we present a technically correct, numerically relevant description of the consequences of rigorous measure theory that is accessible to a wide technical audience. We then describe how we exploit these results to produce a numerical solution technique.

The ingredients of measure theory are: (1) a domain; (2) a σ-algebra; and (3) a (probability) measure. After specifying the domain of interest (Λ), a σ-algebra (ℬΛ) is the class of subsets in the domain for which it makes sense to compute the (probability) measure. In the case of a probability measure, the sets in ℬΛ are called events. The pair (Λ, ℬΛ) defines what is called a measurable space. A measure can only be defined once a domain and a σ-algebra are specified (i.e., once a measurable space is defined). There are typically many possible measures on a given measure space.

The notion of a σ-algebra is often presented as an abstract construction involving infinite sequences of set operations, e.g., intersections, unions, complements, beginning with simple sets. However, this construction is actually grounded in practical numerical approximation. Many natural probability questions lead to the need to compute the probability measure of very complex events. This is carried out by approximating complex events by collections of simpler sets, e.g. generalized rectangles when Λ ⊂ ℝ, and then computing a piecewise constant approximate probability measure with respect to the approximating sets. The central theoretical foundation of measure theory insures this approach works.

Another key ingredient of measure theory is the concept of measurable function, or random variable in the case of a probability measure. In particular, measurable functions such as probability densities are those functions which can be integrated over measurable sets. A measurable function is defined in terms of the behavior of its inverse with respect to the σ−algebras on its domain and range. Specifically, if the preimage (inverse) of any set in the output σ-algebra is a set in the input σ-algebra, then the function is said to be measurable. Examples of measurable functions include the physics-based maps from input parameters to observable data and probability densities which map the measurable spaces into the non-negative real numbers.

Given a measure/probability on a measurable space, a measurable function on this space induces a measure/probability on the output space of the measurable function in terms of the inverses of measurable sets. Thus, measure theory is ideal for both the formulation and solution of the inverse problem for scientific inference using deterministic models even in cases where the inverse of the measurable map is set-valued.

4.2. The generalized contours of an inverse problem

In the most idealized case, a map Q : Λ → 𝒟 assigns distinct output values to distinct input values, so the inverse of the map is defined as a pointwise map from 𝒟 to Λ. However, this ideal case is far too restrictive in general. In particular, physics models often have the property that the domain has higher dimension than the range of the observable outputs defining the map Q, which makes it impossible to assign distinct input values to distinct output values.

In some mathematical fields, such maps are referred to as “non-invertible” and the inverse problem is said to be “ill-posed.” However, the inverse of any map is a well defined concept. The inverse of a map is generally a set-valued function that does not map the range 𝒟 back to Λ in a pointwise sense. Rather the inverse of a map sends 𝒟 to a new space ℒ whose points consist of sets (preimages) in Λ. When the map Q is piecewise differentiable, these set-valued inverses exist as lower dimensional manifolds that we refer to as generalized contours. This is a generalization to higher dimensions of the familiar concept of a contour map for elevations. We show an example in Fig. 2, where we plot some of the contours of the map from ℝ2 to ℝ1,

Figure 2

The left plot shows a map from Λ ⊂ ℝ2 to an interval 𝒟 and corresponding contour curves in Λ corresponding to unique points in 𝒟. The right plot shows the decomposition of Λ into a set of contour curves.

For a map between Λ ⊂ ℝ and 𝒟 ⊂ ℝ with m > d > 1, the generalized contour is a manifold of dimension m − d embedded in Λ. For the problem of determining Manning’s n, a point in the space ℒ simply represents all the possible Manning’s n fields producing the same maximum water height. Moreover, we can decompose Λ as a union of generalized contours corresponding to the points in ℒ (see Fig. 2).

Notions such as continuity and well-posedness for the inverse map must be phrased with respect to points in ℒ, not with respect to points in Λ. For example, the condition or degree of near-ill-posedness of the inverse map is determined by how well the map distinguishes different generalized contours from each other. The concept of condition is explored in more depth in Section 6.

The map Q has a natural restriction Qℒ to ℒ since we can evaluate Qℒ on a generalized contour by evaluating Q at any point on the contour. It follows that Qℒ : ℒ → 𝒟 is a bijection, so is invertible, and it is possible to identify each generalized contour by the unique corresponding value in 𝒟 (see Fig. 3). Contours in a contour map are usually indexed by the altitude. However, it is also possible to find manifolds in Λ that can serve as an indexing set for the generalized contours in ℒ (see Fig. 3 and 4) in the sense that each generalized contour corresponds to a unique point of the indexing set. We call such an indexing set a transverse parametrization (TP)[1]. In the context of the Manning’s n problem, a TP simply identifies a way to distinguish which sets of Manning’s n fields produce distinct water heights.

Figure 3

The left plot shows a standard indexing of the set of generalized contours by the corresponding output values of the map. The middle plot shows the contours indexed by distance along a TP representing ℒ. The profile of the map along the TP is shown in Fig. 4. The right plot shows the nonlinear coordinate system induced by the contours and TP, where xℒ is the distance along the TP and x𝒞 is the distance along the contour corresponding to λ.

Figure 4

The left plot is of the profile of the map along the TP from Fig. 3. The right plot shows the inverse density along the TP corresponding to a uniform density on 𝒟.

We choose any TP to serve as a representation of ℒ, and abusing notation, we call this TP ℒ. With this choice, we obtain a new (nonlinear) “coordinate system” for the original Λ. Any point λ ∈ Λ can be written as λ = [xℒ; x𝒞], where xℒ ∈ ℒ indicates the unique generalized contour 𝒞λ passing through λ and x𝒞 is the unique point in 𝒞λ that distinguishes λ from the other points in 𝒞λ (see Fig. 3).

We note that it is possible to numerically approximate all of these manifolds in the set-valued solution of the inverse problem ([23, 24, 25, 26]). However, this is computationally intensive and it is not necessary for solving the stochastic inverse problem.

4.3. Solving the stochastic inverse problem in the space of generalized contours

If the generalized contours are available, then the solution of the stochastic inverse problem into ℒ is conceptually simple since Qℒ is a bijection. A probability distribution on 𝒟 corresponds to a unique induced probability distribution on ℒ. We let P𝒟 denote the given probability measure on 𝒟, which means that if B is any event in 𝒟, then the probability of B is P𝒟(B). The induced inverse probability measure Pℒ on ℒ is given by for any event A ⊂ Λ. This is the “natural” solution of the stochastic inverse problem since it involves minimal assumptions.

Given a probability density ρ𝒟 on 𝒟, there is a unique induced inverse probability density function ρℒ on ℒ such that The induced probability density on the TP from Fig. 3 corresponding to a uniform probability measure on 𝒟 is shown in Fig. 4.

To interpret the solution of the inverse problem in the space of generalized contours ℒ, we note that an event in ℒ corresponds to a set of generalized contours in Λ, which we call a contour event. We illustrate a contour event in Fig. 5 where a probability measure on ℒ provides the means to compute the probability of sets in Λ comprised of generalized contours. For the Manning’s n problem, this means that we can use a probability density on maximum water height data to determine the probability density on events defined by all possible Manning’s n fields that can generate this data.

Figure 5

The left plot shows the inverse density along the TP corresponding to a uniform density on 𝒟. The probability of the contour event shown in red is the integral of the inverse density over the event on the TP determining the contour event shown in a darker shade. The right plot shows four events in Λ defining the same contour event and hence are assigned the same probability by the density on the TP.

4.4. Solving the stochastic inverse problem in the original domain

The formulation of the parameter domain Λ is an important ingredient in the construction of a physics model. Parameter values correspond to physical conditions that determine model behavior and, generally, different points and events in Λ correspond to different physical conditions even when they result in the same model output. Many important scientific and engineering applications of a model require the ability to explore differences between various sets of parameter values in Λ.

The solution of the stochastic inverse problem in the space of generalized contours provides the means to answer questions in Λ that are formulated in terms of contour events. We can identify arbitrary events A and B in Λ by the smallest contour events containing A and B, i.e., by Q−1(Q(A)) and Q−1(Q(B)), respectively. However, without further assumption, we cannot use the solution of the stochastic inverse problem in ℒ to assign different probabilities to different events A and B in Λ when Q(A) = Q(B) (see Fig. 5).

In the deterministic inverse problem, the model cannot be used to distinguish between different points on the same generalized contour. However, there is a reasonable assumption that allows assigning probabilities to individual events in Λ that are in the same contour event in the case of the stochastic inverse problem. To explain this, we formulate the probability computations in the common form as integrals involving probability densities.

Consider any probability distribution PΛ that is given on Λ in terms of a probability density ρΛ, so that for any event A ⊂ Λ, We introduce notation to exploit the coordinate system induced by the generalized contours. We fix a TP representing ℒ in Λ. For an event A ⊂ Λ, we let π(A) ⊂ ℒ denote the event in ℒ corresponding to the contour event determined by A. (We can compute π(A) as the unique event in ℒ determined by Q−1 (Q(A))). For a point xℒ ∈ ℒ, we let π−1(xℒ) be the generalized contour corresponding to xℒ (see Fig. 6).

Figure 6

The left plot shows the variables used to decompose the integral (4.4). The right plot shows the densities along different contours determined using the standard assumption of a uniform density. Note the decreased height of the density for the generalized contours that are apparently longer.

Changing to the new coordinates, we obtain the iterated integral, where we abuse notation to let x𝒞 denote the coordinates along the generalized contour 𝒞 = π−1(xℒ) corresponding to xℒ. Moreover, the density ρΛ on Λ can be decomposed as where ρℒ is a probability density function on ℒ and ρ𝒞(x𝒞; xℒ) is a conditional probability density conditioned on the generalized contour π−1(xℒ). Substituting, we obtain,

We note that when ρΛ is given, then ρℒ is obtained by restriction of ρΛ to ℒ and ρ𝒞 is obtained by restriction of ρΛ to the generalized contours corresponding to points in ℒ. Vice versa, if we specify ρℒ and ρ𝒞, then ρΛ is determined uniquely.

We emphasize the great theoretical and computational importance of the iterated integral (4.7) for probability measure on Λ. Roughly speaking, (4.7) states the intuitive fact that the probability of an event A in Λ can be computed as the product of the conditional probability of A considered as an event in the probability space of the contour event π(A) and the probability of that contour event. However, (4.5) and (4.6) are not obvious from a mathematical point of view. Technical issues potentially arise because ℒ and the generalized contours are lower dimensional structures in Λ and so have measure zero in the volume measure of Λ[2].

We now return to the stochastic inverse problem. As discussed above, (4.3) determines ρℒ uniquely from the inverse of the model and the imposed probability density ρ𝒟 on 𝒟. However, the model gives information about the conditional probability ρ𝒞 in each generalized contour. In the context of the Manning’s n problem, the probability density of water height data uniquely determines the probability density ρℒ but provides no information about the possible conditional probability ρ𝒞 which would distinguish between different Manning’s n values that give the same height. At this point, we make an assumption.

There is a branch of probability that deals with the assignment of probabilities when there is no information about probabilities of events. In the case of a discrete space, the Principle of Insufficient Reason[3], going back to Bernoulli and Laplace, states that if there is no reason to assume that one outcome is favored over any others, we should assign equal probabilities to all outcomes. The analog in continuous probability is to assign the uniform probability distribution in the situation in which we have no reason to believe that there is nonuniform behavior in the probability.

Thus, we assume that ρ𝒞 is the uniform probability density on each generalized contour With this assumption, (4.7) gives a unique solution of the stochastic inverse problem as a probability density on Λ. The conditional density ρ𝒞 varies with the generalized contour (see Fig. 6), which has a complicated effect on the resulting probability density on Λ (see Fig. 7). In particular, suppose two events A1 and A2 have equal measure in Λ and events π(A1), π(A2) ⊂ ℒ have equal probability Pℒ(π(A1)) = Pℒ(π(A2)) in ℒ. If π(A1) has greater volume than π(A2), then PΛ(A1) < PΛ(A2). This is evident in Fig. 7.

Figure 7

Plots of the inverse density in Λ corresponding to the example shown in Figures 2, 3, 4, and 5 under the standard assumption of uniform probabilities along generalized contours.

We emphasize that any assumed conditional probability density ρ𝒞 determines a unique probability distribution on Λ, and our computational methods work equally well with any assumed density. However, note that any non-uniform density imposes additional geometric structure on the problem that is not determined by the model. For example, assuming a Gaussian distribution on each generalized contour requires specifying a point in each generalized contour as the mean of the Gaussian. The model itself cannot distinguish a specific point in a generalized contour, so we have imposed an additional geometric feature in the problem. For further discussion of the assumption, see [26].

4.5. Comparison to other formulations of stochastic inverse problems

Modern research on the approximate solution of inverse problems is dominated by Bayesian and regularization approaches, and it is impossible to list all the relevant references here. Below, we describe the differences in our approach to some recent work using other approaches in the determination of Manning’s n. For more information, the interested reader should review the cited work and the literature cited therein.

In [41], the optimization approach of [27] is generalized and regularization is used to determine Manning’s n from noisy water height measurements using a shallow water equation model. In the stated inversion algorithm of [41], the inverse problem is formulated as finding a Manning’s n field that minimizes a misfit functional. This requires the specification of a regularization parameter and the objective is to find a specific Manning’s n field for the specified misfit functional. This is a completely different formulation of the inverse problem where the physical map defining the observable data of water heights is used in part to define a separate mapping (the misfit functional) over the parameter space. Regularization approaches such as this have proven useful in identifying a specific input to a model whose corresponding output values are a reasonable match, in some norm, to the observable data. However, the misfit functional, by design, has a completely different contour structure from the physical map we consider in order to produce a so-called “unique” answer to the inverse problem.

In [29, 30], Bayesian approaches are used to determine posterior distributions of Manning’s n from noisy water elevation data. Both of these works discuss methods for using and accelerating MCMC to sample form the posterior distribution and deal with the so-called hyperparameters defining prior distributions. The posterior distribution is defined as a conditional probability density using a likelihood and prior that, under certain assumptions, has an implicit connection to regularization [42]. Furthermore, the likelihood involves the treatment of the physical deterministic map Q as a statistical map where a probability density is implicitly defined over Λ by considering discrepancies in the map Q from observed data. This is fundamentally different than imposing uncertainty in the form of a probability density on the range of Q. The contours of a probability density defined by discrepancies are entirely dependent on the form of the specified density not on the actual physics defining the map Q. Moreover, the posterior is treated typically as an “update” to prior information (e.g., see Fig. 9 in [29] and Figs. 14–16 in [30]). This is often done for the same purpose as regularization approaches, which is to identify a single parameter maximizing the density (corresponding to a minimization of a misfit functional). As with regularization, Bayesian approaches have proven extremely useful at identifying specific input parameters with desired characteristics that explain specific sets of observable output data. However, the contours of the posterior distribution have a completely different interpretation than the contours of the physical map on which we build a non-parametric probability distribution.

We note one consequence of the set-valued nature of the inverse of the forward map is that there are an infinite number of probability structures that will give the same output [23]. Consequently, using the output to quantify the accuracy of a computed inverse solution is a relatively weak metric in comparing these different methods. Further complicating direct comparisons of solutions obtained by these approaches is the fact that the underlying mathematical formulations are fundamentally different leading to solutions with different interpretations.

4.6. Numerical approximation of the inverse probability distribution in Λ

We now adopt (4.8) in the iterated integral (4.7). Formally, this means that we compute PΛ(A) for any event A in Λ by evaluating (4.7). In this section, we review the numerical aspects of this computation and refer to [23, 24, 25, 26] for more details. We finish with an interpretation of the computational algorithm in the context of the Manning’s n problem.

There are two key aspects to approximating an iterated integral over a set, namely dealing with the geometry of the set and evaluating the integrand. We can handle both simultaneously using the standard measure theoretic technique of simple function approximations which are defined as piecewise constant functions on a specified collection of measurable sets.

Generally, evaluating any iterated integral over a set with a non-Cartesian geometry is difficult. Evaluating (4.7) is made more complicated by the generally complex geometries of the TP representing ℒ and the generalized contours, even assuming those structures have been approximated. Moreover, a number of applications involving use of an inverse probability distribution require computing the probability of a large number of events. A way to deal with these issues is to use the fact that an integral of a function f over Λ can be written as a sum of integrals, over a (discretization) partition of Λ consisting of a collection of subdomains that form a nonoverlapping partition of Λ, i.e. Λ = ∪𝒱 and 𝒱 ∩ 𝒱 = ∅ for i ≠ j. We note that each 𝒱 is a measurable set.

With a partition of Λ in hand, we can formally evaluate (4.7) for any event A as a sum of integrals over . We can deal with complicated geometries by approximating (4.7) using either the (outer) sum over the collection of subdomains {𝒱 : A ∩ 𝒱 ≠ ∅} or the (inner) sum over the collection {𝒱 : 𝒱 ⊂ A}, provided that the integrals over the subdomains are easier to evaluate and the differences between A and these collections have small measure.

The second aspect of computing (4.7) is evaluating the integrand, which entails evaluating the inverse density function ρℒ on ℒ defined through the inverse and the uniform probability density function ρ𝒞 along each generalized contour 𝒞. Since these functions are complicated in general, we use piecewise constant approximations defined with respect to partitions of small cells over which the densities are not expected to vary much.

These observations suggest that it might be necessary to construct a number of partitions, namely for Λ, ℒ, and each of the generalized contours. Fortunately, we can construct partitions of Λ that can be used to approximate sets in Λ, ℒ, and each generalized contour event simultaneously. This is based on the fact that any open set in ℝ can be approximated arbitrarily well (in measure) by a collection of generalized rectangles defining “boxes” [40] (see Fig. 8).

Figure 8

Events in Λ are discretized using a Cartesian set of cells. This collection simultaneously approximates events in Λ, ℒ, and the generalized contours 𝒞. We show approximations of an event A, π(A) in ℒ, and the intersection of a contour with A.

Further considerations are required for the numerical evaluation of ρℒ in the part of the integral (4.7) along ℒ and the uniform probability density function ρ𝒞 along each generalized contour 𝒞. Recall that ρℒ is defined as the density induced by the map Q through (4.3), so that evaluating ρℒ formally involves computing integrals ∫ ρ𝒟(y) dy for events A ⊂ ℒ. To do this numerically, we also construct a discretization partition of 𝒟 and compute approximations of the associated probabilities p = ∫ ρ𝒟(y) dy. This yields the approximations such as, We approximate ∫𝒞∩ dx on a given generalized contour 𝒞 by using the ratio of the measure of the cells 𝒱 that intersect 𝒞 ∩ A to the measure of the cells that intersect 𝒞.

Finally, we touch on an important computational point. The approximation of open sets by collections of boxes is theoretically important in any finite dimension, but is computationally important only in low dimensions. The reason is that the number of cells required to partition the unit cell in ℝ increases exponentially with the dimension n, which is a form of the well known “curse of dimensionality”.

The classic way to deal with this dimension dependence is to use a Monte Carlo approximation since the choice of the random sample points is not tied to a Cartesian geometry. However, for the solution of the stochastic inverse problem, it is essential to tie Monte Carlo approximations of integrals to the approximation of events in Λ. This connection is provided in stochastic geometry [43, 44], which studies the properties of stochastic partitions of a domain based on collections of randomly chosen points in Λ. A collection of points in Λ determines a Voronoi tessellation [43] . This is the collection of cells whose sides are, roughly speaking, determined as segments that lie equidistant between neighboring points λ( and λ( (see Fig. 9). A fundamental approximation result is that we can approximate events in Λ using such collections.

Figure 9

We compare approximations of an event A constructed using a Cartesian set of cells and a Voronoi tessellation corresponding to a set of randomly chosen points.

However, constructing Voronoi tesselations in high dimensions is also prohibitively expensive. So the key for practical numerical computations is to construct approximations that exploit the approximation properties of the tesselations implicitly but do not require explicit construction of the cells.

We summarize the numerical approximation for the inverse probabilities in Algorithm 1. With the approximate probabilities {ρ̂Λ,}, we can construct approximate density plots for the inverse density and compute approximations, for example the inner sum, We can also search for regions of highest probability.

Algorithm 1

Numerical Approximation of the Inverse Density

Choose points {λ^(j)}j=1N implicitly defining the partition {𝒱_j}j=1N⊂Λ of Λ.
Assign values Qj = Q(λ(j)) for points λ(j), j = 1, … N.
Choose a discretization partition {I_i}i=1M⊂𝒟 of 𝒟.
Compute approximations pi ≈ ∫Ii ρ𝒟(y) dy, i = 1, …, M.
Let 𝒞i = {j : Qj ∈ Ii}, i = 1, …, M.
Let 𝒪j = {i : Qj ∈ Ii}, j = 1, … N.
Let Vj be an approximate measure of 𝒱j, j = 1, …, N.
for j = 1, …, N do
	Set ρ̂Λ,j = (Vj/∑k∈𝒞𝒪j Vk) p𝒪j.
end

In the context of the Manning’s n problem, the first step of Algorithm 1 corresponds to choosing a set of possible Manning’s n fields. The second step involves solving the model to determine the maximum water heights associated with each of these Manning’s n fields. The third and fourth steps involve computing an approximation to a prescribed probability density on the maximum water height data. The specification of 𝒞 and 𝒪 corresponds to determining the approximating set of all Manning’s n fields that map to the same set of water height data. The for-loop uses these approximations to compute the probability associated to each sample of a Manning’s n field and its implicitly defined Voronoi cell.

5. Defining the parameter domain Λ

We stress the importance of precisely defining the relevant parameter domain Λ. This is generally more complicated than it might first appear, since it involves consideration of physical properties of the system. After presenting an example to illustrate the importance, we discuss the choice of domain for the Manning’s n parameter field.

We first consider the simple model (4.1) plotted in Fig. 2. The model is valid on all of ℝ2, and on ℝ2, the model has the property that the lengths of the contours increases monotonically as the distance from a fixed set 𝒜 containing (.5, .5) and (1.5, 1.5) increases. However, when we restrict the model to Λ = [0, 2] × [0, 2], the effect is to “cut off” parts of the longer contours, and after some critical distance away from 𝒜, the length of the contours contained in Λ actually decreases as the distance from 𝒜 increases. This leads to the increase in the probability density in Λ near the corners (0, 2) and (2, 0) evident in Fig. 7.

This is an example of the complicated effect that the specification of Λ can have on the inverse probability distribution PΛ. In general, defining the appropriate domain Λ is a critical part of formulating the stochastic inverse problem[4].

5.1. A mesoscale representation of Manning’s n

Manning’s n is a field that varies spatially with the land type and physical conditions. High resolution images, see Fig. 10, show the fine scale changes of land type in a typical coastal region.

Figure 10

NOAA C-CAP land coverage (http://www.csc.noaa.gov/crs/lca/gulfcoast.html)[45].

Since the changes in land type are discrete at a sufficiently fine scale, Manning’s n can be represented in terms of coefficients with respect to a basis for a finite dimensional space of functions. Upon substitution of a chosen representation for Manning’s n into a hydrodynamic model (Eq. 2), the coefficients in the representation define the finite dimensional set of parameters input into the model. In this formulation of the forward model (Eq. 2.2) Manning’s n is time-invariant and constant with respect to the flow characteristics within a particular simulation.

However, a pointwise accurate representation of Manning’s n generally requires a very high number of degrees of freedom. In fact, Manning’s n typically varies on a spatial scale that is much finer than the scale of the cells that can be used for numerical discretization of a hydrodynamic model. In other words, a typical discretization cell for a numerical solution of the model contains a variety of land types. Consequently, representations of Manning’s n typically are defined using an upscaling procedure that employs local averaging.

We employ a standard representation. We construct a triangulation of the spatial domain for the model, consisting of a collection of triangles that form a nonoverlapping partition of Ω, i.e. Ω = ∪T and T ∩ T = ∅ for i ≠ j, and represent Manning’s n as a continuous piecewise linear function on Ω that is linear on each triangle T. If we employ the nodal basis ϕ for the continuous piecewise linear functions on the triangulation, we may then define Manning’s n in terms of the set of M nodal values or equivalently as the vector n ∈ ℝ.

It is quite common to encounter computations in which the triangulation used to represent Manning’s n is the same triangulation used to discretize the hydrodynamic model. However, this choice introduces both significant theoretical and computational issues with respect to formulating and solving the hydrodynamic model. Roughly speaking, unless the spatial discretization is sufficiently fine so as to resolve the variation in Manning’s n pointwise, we significantly change the hydrodynamic model being solved with each change of the spatial discretization mesh. One consequence is that it is impossible to carry out convergence studies. We illustrate with a numerical example below in Sec. 8.

Hence, we fix the mesocale triangulation used to represent the Manning’s n field. Moreover, we refine the mesoscale triangulation to construct the discretization triangulation used to compute numerical solutions of the hydrodynamic model. This provides both computational efficiency in evaluating the Manning’s n parameter in solving the hydrodynamic model and predictable numerical behavior in the computed hydrodynamic solutions.

5.2. Defining the parameter domain Λ

The available data for Manning’s n consists of empirically determined intervals of values for common land (bottom type) classifications, e.g., see Table 1. The ranges are assigned to individual cells in a pixelated version of a high resolution image such as Fig. 10. Ranges for values at pixel cells corresponding to underwater locations cannot be observed by such images and may be determined by expert opinion.

Table 1

Example ranges of values of LA-GAP Manning’s n coefficient [46].

Class	Description	Min	Max
23	Water	0.015	0.030
4	Saline Marsh	0.020	0.065
3	Brackish Marsh	0.020	0.070
2	Intermediate Marsh	0.025	0.080
1	Fresh Marsh	0.030	0.085
5	Wetland Forest	0.080	0.160

A typical mesoscale cell in the Manning’s n representation encompasses several pixel cells. To obtain a range for the nodal value of the representation at a mesoscale grid point, a convex average is employed, see [47, 46]. At each node in the mesoscale mesh, we define a rectangle using the maximum and minimum planar coordinates of the centroids of the mesoscale triangles sharing that node, see Fig. 11.

Figure 11

The mesoscale triangles sharing the common node (drawn in red) are displayed with black lines. The corresponding rectangle used for averaging is displayed as a red dashed rectangle. The centroids of the triangles are shown in blue. The pixel cells holding the intervals for the land types are outlined in gray and green.

We compute the nodal Manning’s n value as the average Manning’s n value within the grid scale based on the land classification pixels. Hence, the Manning’s n value at the jth node in the mesh is where a is the number of pixels (outlined in green in Fig. 11) within the grid scale at the jth node for the ith land classification, λ is the Manning’s n value associated with land type i, m is the number of land types, and . This spatial averaging process is a convex linear map from the Manning n values associated with each of the m land classification types to the mesoscale Manning’s n field. Given bounds for each λ, we obtain the domain Λ. We describe a computationally useful parameterization of the mesoscale Manning’s n field over Λ in Sec. 7.

6. Characterizing the condition of the stochastic inverse problem

In this section, we investigate characteristics of the quantities of interest for the stochastic inverse problem that impact the accuracy of numerical solutions. Recall that in the solution of a linear square system of equations, there is a unique solution when the matrix is invertible. However, the accuracy of a numerical solution generally depends on the condition number of the matrix. Analogously, the difficulty in computing accurate numerical solutions of the stochastic inverse problem depends on a “skewness” property of Jacobian of Q that plays the role of a condition number for the stochastic inverse problem.

6.1. Geometrically distinct quantities of interest

We first consider the abstract stochastic inverse problem. Recall that the generalized contours exist as manifolds of dimension m − d in Λ under suitable conditions on the map Q : Λ ⊂ ℝ → 𝒟 ⊂ ℝ. These conditions are discussed in more detail below.

To explain the idea, we first consider an example of a linear map Q represented by a d × m matrix. If Q has full rank d, so the rows of Q are linearly independent, then the generalized contours exist as m − d dimensional hyperplanes in Λ ⊂ ℝ. On the other hand, suppose Q does not have full rank, say the last row is a linear combination of the other rows. Then, each generalized contour is a hyperplane of dimension m − d + 1 and we may as well consider the map Q̃ obtained by deleting the QoI defined by the last row of Q. Thus, we assume that the chosen QoI form a linearly independent set of functionals.

Extending this idea to nonlinear maps is based on local linearization:

Definition 6.1

The component maps of the d dimensional map Q are geometrically distinct (GD) if the Jacobian of Q has full rank at every point in Λ.

This implies that the map obtained by linearizing Q at a point has linearly independent rows. Under this assumption, we can prove that the generalized contours exist as m − d dimensional manifolds and the TP exists as a d dimensional manifolds [26].

6.2. Condition of the numerical solution of the stochastic inverse problem

As with the property of geometrically distinct, the condition of the inverse problem with respect to Q is an issue involving the dependencies among the d vectors comprising Q. In this section, we fix a TP for ℒ and restrict Q to ℒ to get the invertible map Qℒ from ℒ to 𝒟.

We begin by reviewing the relation between determinant and measure (volume). Assume that a set of d-dimensional vectors υ1, ⋯, υ is given. If we consider the linear transformation from ℝ to ℝ defined by the matrix V whose columns are the vectors {υ1, ⋯, υ}, then for any cube A ⊂ ℝ, where V A denotes the image of A under V and μ is the Lebesgue measure (volume) in ℝ.

We relate this to the degree of skewness introduced by the map V. For i ≤ d, let Pa(υ1, υ2, ⋯, υ) denote the parallelepiped determined by m-dimensional vectors υ1, ⋯, υ, We use μ(Pa(υ1, υ2, ⋯, υ)) to denote the measure (volume) of the parallelepiped as an i-dimensional object. Then, det V = μ(Pa(υ1, υ2, ⋯, υ)). Below, we use | · | to denote the standard Euclidean norm. A fundamental decomposition is

Theorem 6.1

[48] Given m-dimensional vectors υ1, ⋯, υ, there exist vectors such that

We note that if , then the image of a generalized cube is a parallelepiped that is significantly skewed in a direction perpendicular to , see Fig. 12.

Figure 12

Illustration of Theorem 6.1.

Returning to (6.1), the map V can decrease volumes both by simple scaling and by inducing skewness. Scaling does not affect the numerical solution of the inverse problem, but skewness does. For this reason, we define a measure of skewness.

Definition 6.2

For vector υ, we define where the definition is independent of the particular representation (6.2). Then we define

By induction, we can find vectors and such that The interpretation is that there exists a generalized rectangle with orthogonal sides of length , |υ| having the same measure as Pa(υ1, υ2, ⋯, υ), see Fig. 12. We can order the vectors so

For simplicity of exposition, we consider the situation in which the quantities of interest almost fails the assumption of geometrically distinct in the sense q1 is almost linearly dependent on q2, ⋯, q. To be precise,

Assumption 1

We assume skew (Qℒ) = skew (Qℒ, q1) ≫ skew (Qℒ, q2) and μ(Pa(q2, ⋯, q))/|q1| = O(1).

Roughly speaking, the image of a generalized cube under V is a parallelepiped with a relatively “fat” base and a small height caused by severe skewness. We set ε = (skew (Qℒ, q1))−1 and V = μ(Pa(q2, ⋯, q)).

Note that the dimension of the sides of the generalized cube cells in a discretization of 𝒟 must be much smaller than to avoid a poor set approximation of 𝒟, see Fig. 13. Thus, we assume that the cells in a discretization of 𝒟 have size γε|q1|, where 0 < γ ≪ 1. The measures of such cells is γε|q1|.

Figure 13

Skewness affects the size of the cells that can be used to discretize 𝒟. An enlargement shows the part of 𝒟 not covered by the set approximation.

Next, we consider the discretization of ℒ. To avoid errors in computation of the inverse images of cells in 𝒟, we assume that each cell in a discretization of ℒ is associated to only one cell in the discretization of 𝒟. Applying the fundamental decomposition to the image of a generalized cube cell of measure σ under Q, we obtain a skew parallelepiped with measure σε |q1| V. Thus, and the measure of the cells in the discretization of ℒ satisfies, This implies that the number of sample values required to compute the inverse distribution is proportional to We conclude that the number of samples required increases dramatically as the skewness of the map Qℒ increases.

In the nonlinear case, we first partition ℒ into a cover of nonoverlapping cells {𝒞} such that on 𝒞, Qℒ(λ) ≈ J (x)λ, for some point x ∈ 𝒞. We then apply the linear result to J (x) to find that a discretization of 𝒞 requires the number of cells to be proportional to The conclusion is that the number of samples implicitly defining the cells partitioning Λ in Algorithm 1 is proportional to This arises entirely from the skewness induced by the map Qℒ.

We illustrate this argument with a simple example. Consider the linear map Q with matrix where 0 < ε ≤ 2, applied on the domain Λ = [0, 1] × [0, 1] × [0, 1]. The map Q is full rank, but it becomes closer to being deficient as ε approaches zero. The range 𝒟 of Q is the parallelepiped in ℝ2 defined by the nodes {(0, 0), (1, 1), (0, ε), (−1,−1 + ε)}, see Fig. 14. As ε decreases, the map Q increasingly skews the image of a cube.

Figure 14

We plot the image of Λ = [0, 1] × [0, 1] × [0, 1] under the map (6.5) for two values of ε.

The inverse problem Qλ = q has solutions consisting of line segments in Λ that are perpendicular to the coordinate plane (λ1, λ2). We can take ℒ = [0, 1]×[0, 1] in the coordinate plane (λ1, λ2). On ℒ, Q reduces to the matrix Evaluating (6.4), we conclude that the number of samples required in Algorithm 1 of the stochastic inverse problem is proportional to γ−2ε−1.

6.3. Condition of the forw ard prediction problem

The measure of the range of the output of the forward prediction model depends on the measure of the domain of possible input values. In general, a range with larger measure corresponds to a large range of possible outcomes, and hence less precision in predictions of the outcome. In the scientific inference problem, the measure of the domain of possible input values, determined by the stochastic inverse problem, also depends on the skewness of the model. If the model for the stochastic inverse problem has large skewness, then det J is small. This implies that is large, for any small cell B in 𝒟, where the Jacobian is evaluated at a point in B. Thus, even if the observed values for the output of the model are confined to a range of small measure, the corresponding set of inverse values has large measure. This in turn implies that the measure of the corresponding range of the predication may be large.

7. Computational issues

We fix the mesoscale representation of Manning’s n field with M nodes in the mesh and assume there are m land classification types. We use n to denote the values of the mesoscale representation at the nodes. We define m linear ℝ maps from a value associated with each land classification type in the domain to the mesoscale Manning’s n field, so that Thus, b defines the relative contribution of the Manning’s n value associated to land classification i at each node in the mesoscale mesh. Eq. (5.1) implies that the jth component of b is equal to a, where a is the number of pixels for the ith land classification type within the grid scale of the jth node, and N is the total number of pixels within the grid scale of the jth node. This implies that defines a parameterization of the mesoscale Manning’s n field n in (5.1). Fig. 15 shows an example of the parameterization vectors for three land types with λ1 = 0.142, λ2 = 0.161, and λ3 = 0.012 for an inlet domain.

Figure 15

Parameterization vectors b1 (top left), b2 (top right) and b3 (bottom left) for a mesoscale representation of the Manning’s n field for an inlet domain. The bottom right plot shows the representation defined by 0.142b1 + 0.161b2 + 0.012b3. Note that the inlet area (right-hand side) is considered a shallow area and two land classification types are heavily mixed whereas the deeper (left-hand side) region assigns a single land classification type with a lower Manning’s n value as is common.

We construct the parameterization vectors with a Python wrapped version of GridData [47]. We developed the Python packages PolyADCIRC [49] and BET [50] to sample Λ efficiently either on regular grids, with uniform random sampling, or adaptively. Specifically, PolyADCIRC runs batches of P(arallel)ADCIRC simulations, where the number of simulations per batch and number of processors is user determined. Load balancing is handled automatically. The BET package then processes the results using Algorithm 1 to compute the approximation to PΛ and provide visualizations of results. See Appendix Appendix A for more details.

8. Numerical experiments

In the following set of experiments, we investigate the selection of “effective” quantities of interest for determining the Manning’s n parameters in a hydrodynamic model based on observations of maximum water elevations from a number of possible observation stations.

8.1. The physical domain

We consider an idealized inlet with sloped bathymetry, see Fig. 16. The left boundary is an open ocean boundary with a M2 tidal amplitude of 1.2 [m2/s] entering the domain normal to the boundary. The remainder of the boundary is a land boundary (u · n = 0 and no slip). An earthen jetty acting as an obstacle to tidal flow extends from y = 1500 [m] along the y-coordinate into the domain between x = 1420 [m] and 1580 [m].

Figure 16

Left: Bathymetry (top plot) of the physical domain Ω (with finite element mesh shown on the bottom plot) in Eq. (2) with observation stations marked by circles. Right: The x, y-coordinates of observation stations in Ω for observing QoI. Observations of a QoI from the ith observation station are denoted q(λ). We examine particular subsets of all the possible QoI maps, e.g., Q(λ) = (q1(λ), q5(λ), q12(λ)) or Q(λ) = (q1(λ), q7(λ)).

8.2. Input parameters

We consider two cases. We first fix the length of the earthen jetty to be y = −1050 [m]. Therefore, Λ = [0.07, 0.15] × [0.1, 0.2] is defined by λ1 and λ2 for the parameterization vectors of the Manning’s n field shown in Fig. 15. In the second case, we also allow the length of the earthen jetty in the y-direction to vary in [−1500, 1500] (m), giving Λ = [−1500, 1500] × [0.07, 0.15]× [0.1, 0.2], where we let λ1 denote the length in the y-direction of the earthen jetty and λ2 and λ3 denote the Manning’s n values for land classification types 1 and 2, respectively.

8.3. The potential quantities of interest

We construct the QoI vector with components given by the maximum elevation recorded at a selection of 2 from 12 possible observation stations shown in Fig. 16. The condition of the resulting QoI depends strongly on the physical locations of the stations in the pair and the purpose of the experiment is to investigate the condition of the QoI for various choices. For simplicity, we always fix the first component of the various QoI to be the maximum elevation q1 at station 1 and allow the second component to vary among stations 2 – 12.

8.4. Resolution of the Manning’s n representation and numerical solution

As mentioned, a mesoscale representation of Manning’s n is often defined on the finite element mesh. Since we are modeling the flows at a scale where both the assumptions justifying the validity of the SWE and the finite element mesh has adequately resolved variability in parameters, significant numerical errors can arise with this choice. In our approach, we use a finite element mesh obtained by refinement of the mesoscale mesh, see Fig. 16. The resulting QoI maps are noticeably smoother than those obtained using the common choice, see Fig. 17. Using the refined finite element mesh results in QoI values shifted by approximately 0.5 [m] throughout the entire parameter domain, and the slopes of the contours are substantially different. The unrefined finite element mesh does not sufficiently resolve the earthen jetty nor the varaiblity in the parameters which is on a similar physical scale as the earthen jetty. Refining the finite element mesh aqdeuately resolves these features.

Figure 17

The map Q(λ) = q1(λ) of the maximum elevation measured at station 1 computed on a finite element mesh defined by the mesh used for the mesoscale representation of Manning’s n (left) and on the refined finite element mesh shown in Fig. 16 (right). Note that scales, smoother response, and change in geometry of the contours resulting from using the refined finite element mesh.

In the first case where Λ = [0.07, 0.15]×[0.1, 0.2], we compute the QoI values on a regular 21 × 21 grid of Λ. For the second case where Λ = [−1500, 1500] × [0.07, 0.15] × [0.1, 0.2], we compute the QoI values on a regular 21 × 21 × 21 grid of Λ. Both sets of these simulations were run on Lonestar Linux Cluster at the Texas Advanced Computering Center (TACC) at The University of Texas at Austin. Each compute node of Lonestar contains two processors for a total of 12 cores and 24GB of memory per node [51]. A PADCIRC simulation of a single sample λ ∈ Λ running on a single node with four MPI tasks takes about 90 seconds to run to completion.

8.5. Results for a fixed jetty length

Inverse problem

In Fig. 18, we show representative plots of 𝒟 ≔ Q(Λ) = (q1(Λ), q(Λ)) for k = 2 and k = 6. Recalling Fig. 14, the left plot shows an example of a badly conditioned QoI map relative to the well-conditioned QoI shown in the right plot. The QoI for k = 1 and 3 are similarly ill-conditioned while the QoI for the rest are more or less well-conditioned.

Figure 18

The estimated output domain 𝒟: = Q(Λ) = (q1(Λ), q(Λ)) for k = 2 (left) and k = 6 (right). Comparing to Fig. 14, we see that the left plot indicates a “bad” condition for the inverse problem and the right plot indicates a “good” condition for the inverse problem. Three reference QoI values associated with (λ1, λ2) = (0.107, 0.106), (0.075, 0.121), and (0.0781, 0.168) are marked on each plot.

We define and solve a total of 6 stochastic inverse problems for different ρ𝒟 determined from the three choices of reference QoI values and the two QoI maps shown in Fig. 18. For each problem, we choose ρ𝒟 to be a piecewise constant function with support on small rectangles centered at the reference QoI values. We solve the inverse problem to obtain p̂Λ, for the Voronoi cells defined by regular sampling of the QoI map. We use kernel density estimation to smooth the resulting density approximation.

Figures 19–21 show the resulting approximate probability densities ρΛ and marginals for each of the stochastic inverse problems. We see that using the better conditioned QoI map results in densities with smaller support in Λ, which can have a large effect on the precision of predictions. We note the similarities in the contours of the densities resulting from inverting ρ𝒟 for Q = (q1, q2) to the contours of the map defined only by q1 seen in the righthand plot of Fig. 17. Since q2 coincides very closely with q1, the resulting inverse densities appear as if the geometric information contained in the generalized contours of only q1 (or q2) exclusively was used to construct the inverse density. In Fig. 22, we show the inverse density and marginals resulting from an inversion using only Q = (q1) where the density ρ𝒟 was defined by the projection of the density used for the map Q = (q1, q2) with reference parameter λ = (0.107, 0.106). Comparing Fig. 22 to Fig. 19, we observe that the addition of QoI q2 fails to significantly change the inverse density obtained by inverting q1.

Figure 19

Plots of ρΛ (left) and marginals ρλ (middle) and ρλ (right) inverting ρ𝒟 using Q = (q1, q2) (top) and Q = (q1, q6) (bottom). Here, ρ𝒟 is defined as a uniform density on a small rectangular box centered at the reference QoI values associated with λ = (0.107, 0.106). The reference value and its components are illustrated by a black circle in the ρΛ plots and vertical lines in the marginal plots.

Figure 21

Plots of ρΛ (left) and marginals ρλ (middle) and ρλ (right) inverting ρ𝒟 using Q = (q1, q2) (top) and Q = (q1, q6) (bottom). Here, ρ𝒟 is defined as a uniform density on a small rectangular box centered at the reference QoI values associated with λ = (0.0781, 0.168). The reference value and its components are illustrated by a black circle in the ρΛ plots and vertical lines in the marginal plots.

Figure 22

Plots of ρΛ (left) and marginals ρλ (middle) and ρλ (right) inverting ρ𝒟 using Q = (q1). Here, ρ𝒟 is defined as a uniform density on a small interval defined by projecting the rectangular box used to ρ𝒟 for Q = (q1, q2) centered at the reference QoI values associated with λ = (0.107, 0.106). The reference value and its components are illustrated by a black circle in the ρΛ plots and vertical lines in the marginal plots.

Predictions

The water entering the inlet flows around the bottom of the earthen jetty located at y = −1050 [m] and extending from x = 1420 [m] to 1580 [m] (see Fig. 16). An interesting and important forecasting problem is to determine the time of inundation of critical physical locations within the domain, e.g. the time of inundation near or on physical barriers. With the same model setup, we consider the goal of predicting the time of inundation of points near the bottom of the earthen jetty. Specifically, we consider predicting the time of inundation at two points located at (x1, y1) = (1593.75, −1087.5) and (x2, y2) = (1593.75, −1012.5) corresponding to the nearest nodal points of the finite element mesh to the right of the earthen jetty and equally spaced below and above the bottom of the jetty.

We consider the predicted times of inundation for the Manning’s n parameters defining the regions in Λ containing the largest density values accounting for 95% of all the probability. We also forecast the time of inundation using the reference parameter value. We summarize the results in Table 2 where the time is written as hours:minutes:seconds referenced to the initial model run time. We see that using the better conditioned QoI map for determining an inverse density results in prediction intervals substantially smaller than the poorly conditioned QoI map. At (x1, y1) and (x2, y2) the prediction intervals from the poorly conditioned QoI map are approximately 126% and 794% larger, respectively, than for the better conditioned QoI map.

Table 2

The intervals of prediction for the time of inundation (shown as hours:minutes:seconds) at locations (x1, y1) = (1593.75, −1087.5) (first row) and (x2, y2) = (1593.75, −1012.5) (second row). The interval of prediction was computed by propagating a 95% probability region determined from the densities in Fig. 19. The intervals of prediction are distinguished by different probability regions determined by inverting a density of Q = (q1, q2) (second column) and Q = (q1, q6) (third column). The last column shows the reference time corresponding to the reference parameter that was used in determining the densities on Q.

(x, y)	prediction interval 1	prediction interval 2	Reference time
(x1, y1)	[16:36:50, 19:08:56]	[18:01:36, 19:08:56]	18:36:58
(x2, y2)	[26:35:38, 27:18:52]	[27:14:24, 27:19:14]	27:17:08

8.6. Results for varying jetty length

We let Q = (q1, q5). As before, we define ρ as a uniform density on a small rectangular box centered at a reference QoI value associated with the reference parameter λ = (−600.0, 0.073, 0.119). The marginal density plots are shown in Fig. 23. All of the QoI maps are poorly conditioned, yielding density plots almost identical to those shown in Fig. 23. For example, using Q = (q1, q5, q2) produces what appears to be a 2-D manifold embedded in the set defined by q1(Λ) × q5(Λ) × q2(Λ). In other words, it appears that q2 can be written as a function of q1 and q5 implying it adds almost no useful information for the inverse problem.

Figure 23

Top: Plots of 2-D marginals ρλ (left), ρλ (middle), and ρλ (right). Bottom: Plots of 1-D marginals ρλ (left), ρλ (middle), and ρλ (right). Here, ρ𝒟 is defined as a uniform density on a small rectangular box centered at the reference QoI values associated with λ = (−600.0, 0.073, 0.119). The reference value and its components are illustrated by a black circle in the top plots and vertical lines in the bottom plots.

Since the QoI map is 2-D while Λ is 3-D, the generalized contours are now 1-D curves embedded in the 3-D parameter space. We observe the affects of the generalized contour geometry on the density plots where any increase in the support of a marginal density is a consequence of the increased dimension for the parameters. The QoI map appears to be most sensitive to the length of the earthen jetty given the shape of the marginal densities (i.e., the normalized tangent vectors of the 1-D generalized contour curves have smallest component in the λ1 direction). This is not surprising since the length of this jetty has the greatest affect on restricting the flow of water that reaches station 5 (and thus the value of q5, see Fig. 16).

9. Conclusions

We formulated and numerically solved a stochastic inverse problem involving spatially heterogeneous Manning’s n fields using maximum water elevation data obtained from the ADCIRC model. A novel measure-theoretic framework and computational algorithm was used based on the author’s previous work [23, 24, 25, 26]. However, this previous work did not address the condition of the inverse problem defined in terms of the skewness of contour events. This issue was explored thoroughly in this work. In Section 6, a numerical analysis demonstrated how poor conditioning leads to more samples in Algorithm 1 in order to accurately estimate events in the parameter space. In the numerical experiments, the condition of the inverse problem was explored in the context of the effect on the scientific inference problem. Specifically, poorly chosen quantities of interest lead to a solution of the inverse problem such that predictions based on this solution have reduced precision.

10. Future work

We have demonstrated the utility and highlighted some of the challenges of solving the stochastic inverse problem for quantifying uncertainty in Manning’s n coefficients with ADCIRC within a measure-theoretic framework. We plan to develop goal-oriented adaptive sampling techniques to obtain reasonable estimates of PΛ as we progress to higher dimensional probability spaces. Since many computational models are also computationally expensive another goal of the adaptive sampling is the produce reasonable estimates of PΛ for specific events of physical importance given limited samples. This will require sampling strategies that utilize computed estimates of skewness to optimally place samples within specified events. We also are investigating algorithms for determining the optimal sets of QoI from large data sets (e.g. as results from time series data) using skewness metrics. We plan to apply these mesure-theoretic parameter estimation techniques to a hurricane simulation using a subdomain implementation of ADCIRC on a complex physical domain. We will focus on coastal areas vulnerable to hurricanes (e.g. Southern Louisiana, Southeastern Texas, or the New York and New Jersey coasts). Hurricane simulations on meshes fine enough to resolve inundation are computationally expensive. Thus, we plan to employ a recently available subdomain implementation of ADCIRC [52, 53, 54] to reduce simulation time and allow us to focus on specific areas of interest rather than on the significantly larger domain required for hurricane simulations.