For intelligent ocean exploration and sustainable ocean utilization, the need for smart autonomous underwater vehicles (AUVs), surface craft, and small aircraft is rapidly increasing. Creating time-optimal navigation routes for these vehicles has wide-ranging applications, including ocean data collection, transportation and distribution of goods, naval operations, search and rescue, detecting marine pollution, ocean cleanup, conservation, and solar-windwave energy harvesting. In this thesis, we employ the Massachusetts Institute of Technology– Multidisciplinary Simulation, Estimation, and Assimilation Systems (MIT-MSEAS) time-optimal and hazard-time-optimal path planning theory and schemes based on exact Hamilton–Jacobi partial differential equations (PDEs) and Level Set methods. We apply this methodology to ocean gliders and floats during several real-time sea experiments—the Mini-Adaptive Sampling Test Run (MASTR) and Grand Adaptive Sampling Experiment (GRASE) in the Gulf of Mexico, and the New England Seamounts Acoustic (NESMA) experiment in the North Atlantic. Using the MIT-MSEAS multi-resolution ocean modeling and data assimilation system to provide deterministic and probabilistic ocean current forecasts, we compute time-reachable sets as well as time-optimal paths for a variety of ocean vehicle missions. The governing differential equations for reachability analysis and time-optimal path planning were numerically integrated in real time, forced by our large-ensemble ocean forecasts. We illustrated deterministic and probabilistic forward reachability analyses, glider recovery planning, time-optimal routing for gliders in distress, and planning of future glider and float deployments. Results show that the actual paths of gliders were contained within our reachable set forecasts and in accord with the dynamic reachability fronts. These forecasts were successfully employed for glider recovery and informed strategic decisions for future missions. Additionally, we demonstrated the ability to incorporate risk such as severe weather or vessel traffic into hazard-time-optimal path planning for simulated collaborative air-sea drone missions. Overall, the integration of data-driven multi-resolution ocean modeling with exact reachability theory and numerical schemes enables principled, operationally relevant path planning for diverse ocean missions.

Forecasting circulation in the Gulf of Mexico requires an explicit treatment of uncertainty associated with the Loop Current and its eddies, whose geometry and timing can fluctuate irregularly and lead to chaotic deterministic forecasts. Building on the dynamically orthogonal (DO) methodology for evolving low-rank stochastic representations and on efficient DO numerical schemes for geophysical fluid flows, this thesis develops and assesses massive probabilistic Primitive Equation (PE) hindcasts for the Gulf using the Dynamically Orthogonal Primitive Equations (DO–PE) framework as implemented for realistic ocean dynamics in previous MIT-MSEAS studies. The workflow extracts a time-dependent stochastic subspace from a balanced MIT MSEAS PE ensemble via singular-value decomposition, represents the initial non-Gaussian coefficient cloud with Gaussian mixture models, and subsequently evolves the DO–PE mean, modes, and coefficients under dynamics, numerics, and forcings consistent with the MIT MSEAS PE modeling system.

A 12-day hindcast simulation experiment spanning 28 May–8 June 2015 quantifies skill and convergence across truncations, with weak-type tests (means, standard deviations, kernel-density marginals) and strong-type tests against matched full-order realizations started from identical initial states. Consistent patterns emerge. Uncertainty concentrates along the Loop Current jet, the Yucatán inflow, and eddy peripheries. For weak convergence, as the retained dynamic modes increase from 15 to 60, standard-deviation maps sharpen and expand coherently along these dynamically active features, and the statistics indicate convergence with the normalized RMSEs for both mean and standard deviation fields decreasing in a largely monotonic fashion. At depth and for sea-surface height, late-time mean-error behavior can become mildly non-monotonic, indicating sensitivity to mode allocation among variables. In strong-convergence experiments, DO–PE reconstructions initialized at coefficient quantiles closely track the corresponding full-order trajectories: pathwise misfits remain modest, organize along shear zones, and their RMSE time series lie below persistence and within the envelopes implied by the weak-type spread, reinforcing that truncation primarily filters small-scale content while preserving trajectory-level evolution over the 10–12-day window.

Together, these results demonstrate a practical, reproducible pipeline for massive probabilistic forecasting in the Gulf of Mexico that respects PE dynamics while quantifying and localizing forecast uncertainty in flow-dependent ways (details, configuration, and figures in Chapters 3–4). This thesis also introduces dynamic web pages for the interactive visualization of DO–PE output, facilitating the inspection of mean fields, modes, and standard deviations over time in Chapter 5.

Data-driven, deep-learning modeling frameworks have been recently developed for forecasting time series data. Such machine learning models may be useful in multiple domains including the atmospheric and oceanic ones, and in general, the larger fluids community. The present work investigates the possible effectiveness of such deep neural operator models for reproducing and predicting classic fluid flows and simulations of realistic ocean dynamics. We first briefly evaluate the capabilities of such deep neural operator models when trained on a simulated two-dimensional fluid flow past a cylinder. We then investigate their application to forecasting ocean surface circulation in the Middle Atlantic Bight and Massachusetts Bay, learning from high-resolution data-assimilative simulations employed for real sea experiments. We confirm that trained deep neural operator models are capable of predicting idealized periodic eddy shedding. For realistic ocean surface flows and our preliminary study, they can predict several of the features and show some skill, providing potential for future research and applications.

Accurate numerical models are essential to predict the complex evolution of rapidly changing sea ice conditions and study impacts on climate and navigation. However, sea ice models contain uncertainties associated with initial conditions and forcing (wind, ocean), as well as with parameter values, functional forms of the constitutive relations, and state variables themselves, all of which limit predictive capabilities. Due to the multiple types and scales of sea ice and the complex nonlinear mechanics and high dimensionality of differential equations, efficient ocean and sea ice probabilistic modeling, Bayesian inversion, and machine learning are challenging. In this work, we implement a deterministic 2D viscoplastic sea ice solver and derive and implement new sea ice probabilistic models based on the dynamically orthogonal (DO) equations.

We focus on the stochastic two-dimensional sea ice momentum equations with nonlinear viscoplastic constitutive law. We first implement and verify a deterministic 2D viscoplastic sea ice solver. Next, we derive the new stochastic Sea Ice Dynamically Orthogonal equations and develop numerical schemes for their solution. These equations and schemes preserve nonlinearities in the underlying spatiotemporal dynamics and evolve the non-Gaussianity of the statistics. We evaluate and illustrate the new stochastic sea ice modeling and schemes using idealized stochastic test cases. We employ two stochastic test cases with different types of sea ice: ice sheets and frozen ice cover with uncertain initial velocities. We showcase the ability to evolve non-Gaussian statistics and capture complex nonlinear dynamics efficiently. We study the convergence to the physical discretization, and stochastic convergence to the stochastic subspace size and coefficient samples. Finally, we assess and show significant computational and memory efficiency compared to the direct Monte Carlo method.

Numerical modelling of ocean physics is essential for multiple applications such as scientific inquiry and climate change but also renewable energy, transport, autonomy, fisheries, water, harvesting, tourism, communication, conservation, planning, and security. However, the wide range of scales and interactions involved in ocean dynamics make numerical modelling challenging and expensive. Many regional ocean models resort to a hydrostatic (HS) approximation that significantly reduces the computational burden. However, a challenge is to capture and study local ocean phenomena involving complex dynamics over a broader range of scales, from regional to small scales, and resolving nonlinear internal waves, subduction, and overturning. Such dynamics require multi-resolution non-hydrostatic (NHS) ocean models. It is known that the main computational cost for NHS models arises from solving a globally coupled elliptic PDE for the NHS pressure. Optimally reducing these costs such that the NHS dynamics are resolved where needed is the motivation for this work.

We propose a new multi-dynamics model to decompose a domain into NHS and HS dynamic regions and solve the corresponding models in their subdomains, reducing the cost associated with the NHS pressure solution step. We extend a high-order NHS solver developed using the hybridizable discontinuous Galerkin (HDG) finite element methodology by taking advantage of the local and global HDG solvers for combining HS with NHS solvers. The multi-dynamics is derived, and the first version is implemented in the HDG framework to quantify computational costs and evaluate accuracy using several analyses. We first showcase results on Rayleigh Taylor instability-driven striations to evaluate computational savings and accuracy compared to the standard NHS HDG and finite-volume solvers. We highlight and discuss sensitivities and performance. Finally, we explore parameters that can be used to identify domain regions exhibiting NHS behaviour, allowing the algorithm to dynamically evolve the NHS and HS subdomains.

An important challenge in the problem of producing accurate forecasts of multiscale dynamics, including but not limited to weather prediction and ocean modeling, is that these dynamical systems are chaotic in nature. A hallmark of chaotic dynamical systems is that they are highly sensitive to small perturbations in the initial conditions and parameter values. As a result, even the best physics-based computational models, often derived from first principles but limited by varied sources of errors, have limited predictive capabilities for both shorter-term state forecasts and for important longer-term global characteristics of the true system. Observational data, however, provide an avenue to increase predictive capabilities by learning the physics missing from lower-fidelity computational models and reducing their various errors. Recent advances in machine learning, and specifically data-driven knowledge-based prediction, have made this a possibility but even state-of-the-art techniques in this area have not been able to produce short-term forecasts beyond a small multiple of the Lyapunov time of the system, even for simple chaotic systems such as the Lorenz 63 model. In this work, we develop a training framework to apply neural ordinary differential equation-based (nODE) closure models to correct errors in the equations of such dynamical systems. We first identify the key training parameters that have an outsize effect on the learning ability of the neural closure models. We then develop a novel learning algorithm, broadly consisting of adaptive tuning of these parameters, designing dynamic multi-loss objective functions, and an error-targeting batching process. We evaluate and showcase our methodology to the chaotic Balance Equations in an array of increasingly difficult learning settings: first, only the coefficient of one missing term in one perturbed equation; second, one entire missing term in on perturbed equation; third, two missing terms in two perturbed equations; and finally the previous but with a perturbation being two orders of magnitude larger than the state, thereby resulting in a completely different attractor. In each of these cases, our new multi-faceted training approach drastically increases both state-of-the-art state predictability (up to 15 Lyapunov times) and attractor-reproducibility. Finally, we validate our results by comparing them with the predictability limit of the chaotic BE system under different magnitudes of perturbations.

Onboard forecasting and data assimilation are challenging but essential for unmanned autonomous ocean platforms. Due to the numerous operational constraints for these platforms, efficient adaptive reduced-order models (ROMs) are needed. In this thesis, we first review existing approaches and then develop a new adaptive Dynamic Mode Decomposition (DMD)-based, data-driven, reduced-order model framework that provides onboard forecasting and data assimilation capabilities for bandwidth-disadvantaged autonomous ocean platforms. We refer to the new adaptive ROM as the incremental, stochastic Low-Rank Dynamic Mode Decomposition (iLRDMD) algorithm. Given a set of high-fidelity and high-dimensional stochastic forecasts computed in remote centers, this framework enables i) efficient and accurate send and receive of the high-fidelity forecasts, ii) incremental update of the onboard reduced-order model, iii) data-driven onboard forecasting, and iv) onboard ROM data assimilation and learning. We analyze the computational costs for the compression, communications, incremental updates, and onboard forecasts. We evaluate the adaptive ROM using a simple 2D flow behind an island, both as a test case to develop the method, and to investigate the parameter sensitivity and algorithmic design choices. We develop the extension of deterministic iLRDMD to stochastic applications with uncertain ocean forecasts. We then demonstrate the adaptive ROM on more complex ocean fields ranging from univariate 2D, univariate 3D, and multivariate 3D fields from multi-resolution, data-assimilative Multidisciplinary Simulation, Estimation, and Assimilation Systems (MSEAS) reanalyses, specifically from the real-time exercises in the Middle Atlantic Bight region. We also highlight our results using the Navy’s Hybrid Coordinate Ocean Model (HYCOM) forecasts in the North Atlantic region. We then apply the adaptive ROM onboard forecasting algorithm to interdisciplinary applications, showcasing adaptive reduced-order forecasts for onboard underwater acoustics computations and forecasts, as well as for exact time-optimal path-planning with autonomous surface vehicles.

For stochastic forecasting and data assimilation onboard the unmanned autonomous ocean platforms, we combine the stochastic ensemble DMD method with the Gaussian Mixture Model – Dynamically Orthogonal equations (GMM-DO) filter. The autonomous platforms can then perform principled Bayesian data assimilation onboard and learn from the limited and gappy ocean observation data and improve onboard estimates. We extend the DMD with the GMM-DO filter further by incorporating incremental DMD algorithms so that the stochastic ensemble DMD model itself is updated with new measurements. To address some of the inefficiencies in the first combination of the stochastic ensemble DMD with the GMM-DO filter, we further introduce the GMM-DMD algorithm. This algorithm not only uses the stochastic ensemble DMD as a computationally efficient forward model, but also employs the existing decomposition to fit the GMM to and perform Bayesian updates on. We demonstrate this incremental stochastic ensemble DMD with GMM-DO and GMMDMD using a real at-sea application in the Middle Atlantic Bight region. We employ a 300 member set of stochastic ensemble forecasts for the “Positioning System for Deep Ocean Navigation – Precision Ocean Interrogation, Navigation, and Timing” (POSYDON-POINT) sea experiment, and highlight the capabilities of reduced data assimilation using simulated twin experiments.

For intelligent ocean exploration and sustainable ocean utilization, the need for smart autonomous underwater vehicles (AUVs), surface craft, and small aircraft is rapidly increasing. The challenge of creating time-optimal navigation routes for these vehicles has many applications, including ocean data collection, transportation and distribution of goods, naval operations, search and rescue, detecting marine pollution, ocean cleanup, conservation, and solar-wind-wave energy harvesting, among others. In this thesis, we employ the Massachusetts Institute of Technology – Multidisciplinary Simulation, Estimation, and Assimilation Systems (MIT-MSEAS) time-optimal path planning theory and schemes based on exact Hamilton–Jacobi partial differential equation (PDE) and Level Set methods to predict and study the sensitivity of reachable sets and time-optimal trajectories in the Portugal–Azores–Madeira region of the Northern Atlantic, for several types of missions and autonomous ocean vehicles. Specifically, using the MIT-MSEAS multi-resolution ocean modeling and data assimilation system to provide four-dimensional ocean currents in the region, we compute time-reachable sets and time-optimal paths for several missions, and examine the sensitivity to variations in vehicle type, speed, start time, voyage direction, and operating depths. Our real-data-driven multi-resolution simulation study illustrates how navigational paths vary with these parameters, and how ocean dynamics and variability in the Portuguese ocean regions affect the time optimization, as compared to direct voyages in the absence of any ocean currents. We also highlight effects of the Azores and Madeira archipelagos, differences between surface and bottom path planning, interception routes between vehicles of different speeds, and the utilization of arrival time fields in planning. Results showcase how principled path planning, integrating data-driven multi-resolution ocean modeling with exact reachability theory and numerical schemes, can assess the capabilities of ocean vehicles in the Portugal–Azores–Madeira ocean region, by predicting the fastest travel time, expected range, and optimal headings, for varied types of ocean missions.

There are many significant challenges for unmanned autonomous platforms at sea including predicting the likely scenarios for the ocean environment, quantifying regional uncertainties, and updating forecasts of the evolving dynamics using their observations. Due to the operational constraints such as onboard power, memory, bandwidth, and space limitations, efficient adaptive reduced order models (ROMs) are needed for onboard predictions. In the first part, several reduced order modeling schemes for regional ocean forecasting onboard autonomous platforms at sea are described, investigated, and evaluated. We find that Dynamic Mode Decomposition (DMD), a data-driven dimensionality reduction algorithm, can be used for accurate predictions for short periods in ocean environments. We evaluate DMD methods for ocean PE simulations by comparing and testing several schemes including domain splitting, adjusting training size, and utilizing 3D inputs. Three new approaches that combine uncertainty with DMD are also investigated and found to produce practical and accurate results, especially if we employ either an ensemble of DMD forecasts or the DMD of an ensemble of forecasts. We also demonstrate some results from projecting/compressing high-fidelity forecasts using schemes such as POD projection and K-SVD for sparse representation due to showing promise for distributing forecasts efficiently to remote vehicles. In the second part, we combine DMD methods with the GMM-DO filter to produce DMD forecasts with Bayesian data assimilation that can quickly and efficiently be computed onboard an autonomous platform. We compare the accuracy of our results to traditional DMD forecasts and DMD with Ensemble Kalman Filter (EnKF) forecast results and show that in Root Mean Square Error (RMSE) sense as well as error field sense, that the DMD with GMM-DO errors are smaller and the errors grow slower in time than the other mentioned schemes. We also showcase the DMD of the ensemble method with GMM-DO. We conclude that due to its accurate and computationally efficient results, it could be readily applied onboard autonomous platforms. Overall, our contributions developed and integrated stochastic DMD forecasts and efficient Bayesian GMM-DO updates of the DMD state and parameters, learning from the limited gappy observation data sets.

The past decade has seen an increasing use of autonomous vehicles (propelled AUVs, ocean gliders, solar-vehicles, etc.) in marine applications. For the operation of these vehicles, efficient methods for path planning are critical. Path planning, in the most general sense, corresponds to a set of rules to be provided to an autonomous robot for navigating from one configuration to another in some optimal fashion. Increasingly, having autonomous vehicles that optimally collect/harvest external fields from highly dynamic environments has grown in relevance. Autonomously maximizing the harvest in minimum time is our present path planning objective. Such optimization has numerous impactful applications. For instance, in the case of energy optimal path planning where long endurance and low power are crucial, it is important to be able to optimally harvest energy (solar, wind, wave, thermal, etc.) along the way and/or leverage the environment (winds, currents, etc.) to reduce energy expenditure. Similarly, autonomous marine cleanup or collection vehicles, tasked with harvesting plastic waste, oil spills, or seaweed fields, need to be able to plan paths that maximize the amount of material harvested in order to optimize the cleanup or collection process. In this work, we develop an exact partial differential equation-based methodology that predicts harvest-time optimal paths for autonomous vehicles navigating in dynamic environments. The governing differential equations solve the multi-objective optimization problem of navigating a vehicle autonomously in a highly dynamic flow field to any destination with the goal of minimizing travel time while also maximizing the amount harvested by the vehicle. Using Hamilton-Jacobi theory for reachability, our methodology computes the exact set of Pareto optimal solutions to the multi-objective path planning problem. This is completed by numerically solving a reachability problem for the autonomous vehicle in an augmented state space consisting of the vehicle’s position in physical space as well as its harvest state. Our approach is applicable to path planning in various environments, however we primarily present examples of navigating in dynamic ocean flows. The following cases, in particular, are studied. First, we validate our methodology using a benchmark case of planning paths through a region with a harvesting field present in a halfspace, as this case admits a semi-analytical solution that we compare to the results of our method. We next consider a more complex unsteady environment as we solve for harvest-time optimal missions in a quasi-geostrophic double-gyre ocean flow field. Following this, we provide harvest-time optimal paths to the highly relevant issue of collecting harmful algae blooms. Our final case considers an application to next generation offshore aquaculture technologies. In particular, we consider in this case path planning of an offshore moving fish farm that accounts for optimizing fish growth. Overall, we find that our exact planning equations and efficient schemes are promising to address several pressing challenges for our planet.

The present work focuses on applying the Dynamically Orthogonal Primitive Equations (DO-PE) for realistic high-resolution stochastic ocean forecasting in regions with complex ocean dynamics. In the first part, we identify and test a streamlined process to create multi-region initial conditions for the DO-PE framework, starting from temporally and spatially sparse historical data. The process presented allows us to start from a relatively small but relevant set of measured temperature and salinity historical vertical profiles (on the order of hundreds) and to generate a massive set of initial conditions (on the order of millions) in a stochastic subspace, while still ensuring that the initial statistics respect the physical processes, modeled complex dynamics, and uncertain initial conditions of the examined domain. To illustrate the methodology, two practical examples—one in the Gulf of Mexico and another in the Alboran Sea—are provided, along with a review of the ocean dynamics for each region. In the second part, we present a case study of three massive stochastic DO-PE forecasts, corresponding to ensembles of one million members, in the Gulf of Mexico region. We examine the effect of adding more dynamic DO modes (i.e., stochastic dimensions) and show that it tends to statistical convergence along with an enhancement of the uncertainty captured by the DO forecast realizations, both by increasing the variance of already existing features as well as by adding new uncertain features. We also use this case study to validate the DO-PE methodology for realistic high-resolution probabilistic ocean forecasting. We show good accuracy against equivalent deterministic simulations, starting from the same initial conditions and simulated with the same assumptions, setup, and original ocean model equations. Importantly, by comparing the reduced-order realizations against their deterministic counterparts, we show that the errors due to the DO subspace truncation are much smaller and growing slower than the fields themselves are evolving in time, both in the Root Mean Square Error (RMSE) sense as well as in the 3D multivariate ocean field sense. Based on these observations, we conclude that the DO-PE realizations closely match their full-order equivalents, thus enabling massive forecast ensembles with practically low numerical errors at a tractable computational cost.

Though computing power continues to grow quickly, our appetite to solve larger and larger problems grows just as fast. As a consequence, reduced-order modeling has become an essential technique in the computational scientist’s toolbox. By reducing the dimensionality of a system, we are able to obtain approximate solutions to otherwise intractable problems. And because the methodology we develop is sufficiently general, we may agnostically apply it to a plethora of problems, whether the high dimensionality arises due to the sheer size of the computational domain, the fine resolution we require, or stochasticity of the dynamics. In this thesis, we develop time integration schemes, called retractions, to efficiently evolve the dynamics of a system’s low-rank approximation. Through the study of differential geometry, we are able to analyze the error incurred at each time step. A novel, explicit, computationally inexpensive set of algorithms, which we call perturbative retractions, are proposed that converge to an ideal retraction that projects exactly to the manifold of fixed-rank matrices. Furthermore, each perturbative retraction itself exhibits high-order convergence to the best low-rank approximation of the full-rank solution. We show that these high-order retractions significantly reduce the numerical error incurred over time when compared to a naive Euler forward retraction. Through test cases, we demonstrate their efficacy in the cases of matrix addition, real-time data compression, and deterministic and stochastic differential equations.

We develop an exact partial differential equation-based methodology that predicts time-energy optimal paths for autonomous vehicles navigating in dynamic environments. The differential equations solve the multi-objective optimization problem of navigating a vehicle autonomously in a dynamic flow field to any destination with the goal of minimizing travel time and energy use. Based on Hamilton-Jacobi theory for reachability and the level set method, the methodology computes the exact Pareto optimal solutions to the multi-objective path planning problem, numerically solving the equations governing time-energy reachability fronts and optimal paths. Our approach is applicable to path planning in various scenarios, however we primarily present examples of navigating in dynamic marine environments. First, we validate the methodology through a benchmark case of crossing a steady front (a highway flow) for which we compare our results to semi-analytical optimal path solutions. We then consider more complex unsteady environments and solve for time-energy optimal missions in a quasi-geostrophic double-gyre ocean flow field.

Developing accurate and computationally efficient models for ocean acoustics is inherently challenging due to several factors including the complex physical processes and the need to provide results on a large range of scales. Furthermore, the ocean itself is an inherently dynamic environment within the multiple scales. Even if we could measure the exact properties at a specific instant, the ocean will continue to change in the smallest temporal scales, ever increasing the uncertainty in the ocean prediction. In this work, we explore ocean acoustic prediction from the basics of the wave equation and its derivation. We then explain the deterministic implementations of the Parabolic Equation, Ray Theory, and Level Sets methods for ocean acoustic computation. We investigate methods for evolving stochastic fields using direct Monte Carlo, Empirical Orthogonal Functions, and adaptive Dynamically Orthogonal (DO) differential equations. As we evaluate the potential of Reduced-Order Models for stochastic ocean acoustics prediction, for the first time, we derive and implement the stochastic DO differential equations for Ray Tracing (DO-Ray), starting from the differential equations of Ray theory. With a stochastic DO-Ray implementation, we can start from non-Gaussian environmental uncertainties and compute the stochastic acoustic ray fields in a reduced order fashion, all while preserving the complex statistics of the ocean environment and the nonlinear relations with stochastic ray tracing. We outline a deterministic Ray-Tracing model, validate our implementation, and perform Monte Carlo stochastic computation as a basis for comparison. We then present the stochastic DO-Ray methodology with detailed derivations. We develop varied algorithms and discuss implementation challenges and solutions, using again direct Monte Carlo for comparison. We apply the stochastic DO-Ray methodology to three idealized cases of stochastic sound-speed profiles (SSPs): constant-gradients, uncertain deep-sound channel, and a varied sonic layer depth. Through this implementation with non-Gaussian examples, we observe the ability to represent the stochastic ray trace field in a reduced order fashion.

We address the problem of finding the closest matrix Ũ to a given U under the constraint that a prescribed second-moment matrix P̃ must be matched, i.e. Ũ^TŨ=P̃. We obtain a closed-form formula for the unique global optimizer Ũ for the full-rank case, which is related to U by an SPD (symmetric positive definite) linear transform. This result is generalized to rank-deficient cases as well as to infinite dimensions. We highlight the geometric intuition behind the theory and study the problem’s rich connections to minimum congruence transform, generalized polar decomposition, optimal transport, and rank-deficient data assimilation. In the special case of P̃=I, minimum-correction second-moment matching reduces to the well-studied optimal orthonormalization problem. We investigate the general strategies for numerically computing the optimizer, analyze existing polar decomposition and matrix square root algorithms. More importantly, we modify and stabilize two Newton iterations previously deemed unstable for computing the matrix square root, which can now be used to efficiently compute both the orthogonal polar factor and the SPD square root. We then verify the higher performance of the various new algorithms using benchmark cases with randomly generated matrices. Lastly, we complete two applications for the stochastic Lorenz-96 dynamical system in a chaotic regime. In reduced subspace tracking using dynamically orthogonal equations, we maintain the numerical orthonormality and continuity of time-varying base vectors. In ensemble square root filtering for data assimilation, the prior samples are transformed into posterior ones by matching the covariance given by the Kalman update while also minimizing the corrections to the prior samples.

Grand challenges in ocean acoustic propagation are to accurately capture the dynamic environmental uncertainties and to predict the evolving probability density function of stochastic acoustic waves. This is due to the complex ocean physics and acoustics dynamics, nonlinearities, multiple scales, and large dimensions. There are several sources of uncertainty including: the initial and boundary conditions of the ocean physics and acoustic dynamics, the bathymetry and seabed fields; the models parameters; and, the models themselves. In the present work, we start addressing these challenges by deriving, implementing and verifying new optimally-reduced Dynamically Orthogonal (DO) differential equations that govern the propagation of stochastic acoustic waves for underwater sound propagation in an uncertain ocean environment. The developed methodology allows modeling environmental uncertainties in a rigorous probabilistic framework and predicting the uncertainties of acoustic fields, fully respecting the nonlinear governing equations and non-Gaussian statistics of the sound speed and acoustic state variables. The methodology is applied to the standard narrow-angle parabolic equation and is utilized to predict acoustic field uncertainties for three new stochastic idealized test cases: (1) an uncertain Pekeris waveguide with penetrable bottom, (2) an uncertain horizontal interface problem, and (3) an uncertain range-dependent sloping interface problem. For the first case, the solutions of the DO acoustic equations are validated against those obtained using standard Monte Carlo sampling. The second test case showcases results for predicting acoustic field probabilities due to uncertainties in the location of a sound speed channel. For the third test case, the advantages of the DO acoustic equations in predicting uncertainties in complex range-dependent environments are highlighted.

This work focuses on developing efficient and robust implementation methods for hybridizable discontinuous Galerkin (HDG) schemes for fluid and ocean dynamics. In the first part, we compare choices in weak formulations and their numerical consequences. We address details in making the leap from the mathematical formulation to the implementation, including the different spaces and mappings, discretization of the integral operators, boundary conditions, and assembly of the linear systems. We provide a flexible mapping procedure amenable to both quadrature-free and quadrature-based discretizations, and compare the accuracy of the two on different problem geometries. We verify the quadrature-free approach, demonstrating that optimal orders of convergence can be obtained, even on non-affine and curvilinear geometries. The second part of the work investigates the scalability of HDG schemes, identifying memory and time-to-solution bottlenecks. The form of the quadrature-free integral operators is exploited to develop a novel and efficient matrix-free approach to solving the global linear system that arises from HDG discretizations. Additional manipulations to improve numerical robustness are discussed. To mitigate the complexity of the implementation, we provide an automated and computationally efficient verification procedure for the HDG methodologies discussed, using a hierarchical approach to provide diagnostic information and isolate problems. Finally, challenges related to the effective visualization of high-order, discontinuous HDG-FEM data for fluid and ocean applications are illustrated and strategies are provided to address them.

Comprehensive spatiotemporal modeling and forecasting systems for ocean dynamics necessitate robust and efficient data delivery and visualization techniques. The multidisciplinary simulation, estimation, and assimilation systems group at MIT (MSEAS) focuses on capturing and predicting diverse ocean dynamics, including physics, acoustics, and biology on varied scales, thereby developing new methods for multi-resolution ocean prediction and analysis, including data generation and assimilation. The group has primarily used non-interactive ocean plots to visualize its simulated and measured data. Although these maps and sections allow for analysis of ocean physics and the underlying numerical schemes, more interactive maps provide more user control over depicted data, allowing easier study and pattern identification on multiple scales. Integrating static and geospatial data in dynamic visualization creates a heightened viewpoint for analysis, enhances ocean monitoring and prediction, and contributes to building scientific knowledge. This thesis focuses on explaining the motivation behind and the methodologies applied in designing these interactive maps.

This thesis demonstrates the use of exact equations to predict time-optimal mission plans for a marine vehicle that visits a number of locations in a given dynamic ocean current field. The missions demonstrated begin and end in the same location and visit a finite number of locations or waypoints in the minimal time; this problem bears close resemblance to that of the classic “traveling salesman,” albeit with the added complexity of a continuously changing flow field. The paths, or “legs,” between all goal waypoints are generated by numerically solving exact time-optimal path planning level-set differential equations. The equations grow a reachability front from the starting location in all directions. Whenever the front reaches a waypoint, a new reachability front is immediately started from that location. This process continues until one set of reachability fronts has reached all goal waypoints and has returned to the original location. The time-optimal path for the entire mission is then obtained by trajectory backtracking, going through the optimal set of reachability fields in reverse order. Due to the spatial and temporal dynamics, a varying start time results in different paths and durations for each leg and requires all permutations of travel to be calculated. Even though the method is very efficient and the optimal path can be computed serially in real-time for common naval operations, for additional computational speed, a high-performance computing cluster was used to solve the level set calculations in parallel. This method is first applied to several hypothetical missions. The method and distributed computational solver are then validated for naval applications using an operational multi-resolution ocean modeling system of real-world current fields for the complex Philippines Archipelago region. Because the method calculates the global optimum, it serves two purposes. It can be used in its present form to plan multi-waypoint missions offline in conjunction with a predictive ocean current modeling system, or it can be used as a litmus test for approximate future solutions to the traveling salesman problem in dynamic flow fields.

Ocean currents transport a variety of natural (e.g. water masses, phytoplankton, zooplankton, sediments, etc.) and man-made materials (e.g. pollutants, floating debris, particulate matter, etc.). Understanding such uncertain Lagrangian transport is imperative for reducing environmental damage due to natural hazards and for allowing rigorous risk analysis and effective search and rescue. While secondary variables and trajectories have classically been used for the analyses of such transports, Lagrangian Coherent Structures (LCSs) provide a robust and objective description of the important material lines. To ensure accurate and useful Lagrangian hazard scenario predictions and prevention, the first goal of this thesis is to obtain accurate probabilistic prediction of the underlying stochastic velocity fields using the Dynamically Orthogonal (DO) approach. The second goal is to merge data from both Eulerian and Lagrangian observations with predictions such that the whole information content of observations is utilized.

In the first part of this thesis, we develop high-order numerical schemes for the DO equations that ensure efficiency, accuracy, stability, and consistency between the Monte Carlo (MC) and DO solutions. We discuss the numerical challenges in applying the DO equations to the unsteady stochastic Navier-Stokes equations. In order to maintain consistent evaluation of advection terms, we utilize linear centered advection schemes with fully explicit and linear Shapiro filters. We then discuss how to combine the semi-implicit projection method with new high order implicit-explicit (IMEX) linear multi-step and multistage IMEX-RK time marching schemes for the coupled DO equations to ensure further stability and accuracy. We also review efficient numerical re-orthonormalization strategies during time marching. We showcase our results with stochastic test cases of stochastic passive tracer advection in a deterministic swirl flow, stochastic flow past a cylinder, and stochastic lid-driven cavity flow. We show that our schemes improve the consistency between reconstructed DO realizations and the corresponding MC realizations, and that we achieve the expected order of accuracy.

In the second part of the work, we first undertake a study of different Lagrangian instruments and outline how the DO methodology can be applied to obtain Lagrangian variables of stochastic flow maps and LCS in uncertain flows. We then review existing methods for Bayesian Lagrangian data assimilation (DA). Disadvantages of earlier methods include the use of approximate measurement models to directly link Lagrangian variables with Eulerian variables, the challenges in respecting the Lagrangian nature of variables, and the assumptions of linearity or of Gaussian statistics during prediction or assimilation. To overcome these, we discuss how the Gaussian Mixture Model (GMM) DO Filter can be extended to fully coupled Eulerian-Lagrangian data assimilation. We define an augmented state vector of the Eulerian and Lagrangian state variables that directly exploits the full mutual information and complete the Bayesian DA in the joint Eulerian-Lagrangian stochastic subspace. Results of such coupled Eulerian-Lagrangian DA are discussed using test cases based on a double gyre flow with random frequency.

In this work novel high-order hybridizable discontinuous Galerkin (HDG) projection methods are further developed for ocean dynamics and geophysical fluid predictions. We investigate the effects of the HDG stabilization parameter for both the momentum equation as well as tracer diffusion. We also make a correction to our singularity treatment algorithm for nailing down a numerically consistent and unique solution to the pressure Poisson equation with homogeneous Neumann boundary conditions everywhere along the boundary. Extensive numerical results using physically realistic ocean flows are presented to verify the HDG projection methods, including the formation of internal wave beams over a shallow but abrupt seamount, the generation of internal solitary waves from stratified oscillatory flow over steep topography, and the circulation of bottom gravity currents down a slope. Additionally, we investigate the implementation of open boundary conditions for finite element methods and present results in the context of our ocean simulations. Through this work we present the hybridizable discontinuous Galerkin projection methods as a viable and competitive alternative for large-scale, realistic ocean modeling.

Autonomous underwater vehicles (AUVs) are a valuable resource in several oceanic applications such as security, surveillance and data collection for ocean prediction. These vehicles typically travel at speeds comparable to ocean currents, and their movement is significantly affected by these dynamic currents. Further, the speed of currents may vary greatly with depth. Hence, path planning to generate safe and fast vehicle trajectories in such a three-dimensional environment becomes crucial for the successful operation of these vehicles. In addition, many marine vehicles can only move in specific directions and with a speed that is dependent on the direction of travel. Such constraints must be respected in order to plan safe and optimal paths.

Thus, our motivation in this thesis is to study path planning for vehicles with and without motion constraints in three-dimensional dynamic flow-fields. We utilize the time-optimal path planning methodology given by Lolla et al. (2012) for this purpose.

In this thesis, we first review some existing path planning methods (both in two and three-dimensional settings). Then, we discuss the theoretical basis of the rigorous partial differential equation based methodology that is utilized in order to plan safe and optimal paths. This is followed by an elaborate discussion about the application of this methodology to the various types of marine vehicles. We then look at the robust and accurate numerical methods developed in order to solve the governing equations for the path planning methodology with high accuracy in real ocean domains. We illustrate the working and capabilities of our path planning algorithm by means of a number of applications. First we study some benchmark examples with known analytical solutions. Second, we look at more complex flow-fields that analytically model different oceanic flows. Finally, we look at the path planning for different types of marine vehicles in a realistic ocean domain to illustrate the capabilities of the path planning methodology and the developed numerical framework.

This work focuses on developing theory and methodologies for the analysis of material transport in stochastic fluid flows. In a first part, two dominant classes of techniques for extracting Lagrangian Coherent Structures are reviewed and compared and some improvements are suggested for their pragmatic applications on realistic high-dimensional deterministic ocean velocity fields. In the stochastic case, estimating the uncertain Lagrangian motion can require to evaluate an ensemble of realizations of the flow map associated with a random velocity flow field, or equivalently realizations of the solution of a related transport partial differential equation. The Dynamically Orthogonal (DO) approximation is applied as an efficient model order reduction technique to solve this stochastic advection equation. With the goal of developing new rigorous reduced-order advection schemes, the second part of this work investigates the mathematical foundations of the method. Riemannian geometry providing an appropriate setting, a framework free of tensor notations is used to analyze the embedded geometry of three popular matrix manifolds, namely the fixed rank manifold, the Stiefel manifold and the isospectral manifold. Their extrinsic curvatures are characterized and computed through the study of the Weingarten map. As a spectacular by-product, explicit formulas are found for the differential of the truncated Singular Value Decomposition, of the Polar Decomposition, and of the eigenspaces of a time dependent symmetric matrix. Convergent gradient flows that achieve related algebraic operations are provided. A generalization of this framework to the non-Euclidean case is provided, allowing to derive analogous formulas and dynamical systems for tracking the eigenspaces of non-symmetric matrices. In the geometric setting, the DO approximation is a particular case of projected dynamical systems, that applies instantaneously the SVD truncation to optimally constrain the rank of the reduced solution. It is obtained that the error committed by the DO approximation is controlled under the minimal geometric condition that the
original solution stays close to the low-rank manifold. The last part of the work focuses on the practical implementation of the DO methodology for the stochastic advection equation. Fully linear, explicit central schemes are selected to ensure stability, accuracy and efficiency of the method. Riemannian matrix optimization is applied for the dynamic evaluation of the dominant SVD of a given matrix and is integrated to the DO time-stepping. Finally the technique is illustrated numerically on the uncertainty quantification of the Lagrangian motion of two bi-dimensional benchmark flows.

The primary contributions of this thesis include the first stages of development of a 2D, finitevolume, non-hydrostatic, sigma-coordinate code and beginning to apply the Dynamically Orthogonal field equations to study the sensitivity of internal tides to perturbations in the density field. First, we ensure that the 2D Finite Volume (2DFV) code that we use can accurately capture non-hydrostatic internal tides since these dynamics have not yet been carefully evaluated for accuracy in this framework. We find that, for low-aspect ratio topographies, the -coordinate mesh in the 2DFV code produces numerical artifacts near the bathymetry. To ameliorate these staircasing effects, and to develop the framework towards a moving mesh with free-surface dynamics, we have begun to implement a non-hydrostatic sigma-coordinate framework which significantly improves the representation of the internal tides for low-aspect ratio topographies. Finally we investigate the applicability of stochastic density perturbations in an internal tide field. We utilize the Dynamically Orthogonal field equations for this investigation because they achieve substantial model order reduction over ensemble Monte-Carlo methods.

The focus of this research is to study the uncertainties forecast by multi-resolution ocean models and quantify how those uncertainties affect the pressure fields estimated by coupled ocean models. The quantified uncertainty can then be used to provide enhanced sonar performance predictions for tactical decision aides. High fidelity robust modeling of the oceans can resolve various scale processes from tidal shifts to mesoscale phenomena. These ocean models can be coupled with acoustic models that account for variations in the ocean environment and complex bathymetry to yield accurate acoustic field representations that are both range and time independent. Utilizing the MIT Multidisciplinary Environmental Assimilation System (MSEAS) implicit two-way nested primitive-equation ocean model and Error Subspace Statistical Estimation scheme (ESSE), coupled with three-dimensional-inspace (3D) parabolic equation acoustic models, we conduct a study to understand and determine the effects of ocean state uncertainty on the acoustic transmission loss. The region of study is focused on the ocean waters surrounding Taiwan in the East China Sea. This region contains complex ocean dynamics and topography along the critical shelf-break region where the ocean acoustic interaction is driven by several uncertainties. The resulting ocean acoustic uncertainty is modeled and analyzed to quantify sonar performance and uncertainty characteristics with respect to submarine counter detection. Utilizing cluster based data analysis techniques, the relationship between the resulting acoustic field and the uncertainty in the ocean model can be characterized. Furthermore, the dynamic transitioning between the clustered acoustic states can be modeled as Markov processes. This analysis can be used to enhance not only submarine counter detection aides, but it may also be used for several applications to enhance understanding of the capabilities and behavior of uncertainties of acoustic systems operating in the complex ocean environment.

The pressure-correction projection method for the incompressible Navier-Stokes equation is approached as a preconditioned Richardson iterative method for the pressure- Schur complement equation. Typical pressure correction methods perform only one iteration and suffer from a splitting error that results in a spurious numerical boundary layer, and a limited order of convergence in time. We investigate the benefit of performing more than one iteration. We show that that not only performing more iterations attenuates the effects of the splitting error, but also that it can be more computationally efficient than reducing the time step, for the same level of accuracy. We also devise a stopping criterion that helps achieve a desired order of temporal convergence, and implement our method with multi-stage and multi-step time integration schemes. In order to further reduce the computational cost of our iterative method, we combine it with an Aitken acceleration scheme. Our theoretical results are validated and illustrated by numerical test cases for the Stokes and Navier-Stokes equations, using Implicit-Explicit Backwards Difference Formula and Runge-Kutta time integration solvers. The test cases comprises a now classical manufactured solution in the projection method literature and a modified version of a more recently proposed manufactured solution.

The growing use of autonomous underwater vehicles and underwater gliders for a variety of applications gives rise to new requirements in the operation of these vehicles. One such important requirement is optimization of energy required for undertaking missions that will enable longer endurance and lower operational costs. Our goal in this thesis is to develop a computationally efficient, and rigorous methodology that can predict energy-optimal paths from among all time-optimal paths to complete an underwater mission. For this, we develop rigorous a new stochastic Dynamically Orthogonal Level Set optimization methodology. In our thesis, after a review of existing path planning methodologies with a focus on energy optimality, we present the background of time-optimal path planning using the level set method. We then lay out the questions that inspired the present thesis, provide the goal of the current work and explain an extension of the time-optimal path planning methodology to the time-optimal path planning in the case of variable nominal engine thrust. We then proceed to state the problem statement formally. Thereafter, we develop the new methodology for solving the optimization problem through stochastic optimization and derive new Dynamically Orthogonal Level Set Field equations. We then carefully present different approaches to handle the non-polynomial non-linearity in the stochastic Level Set Hamilton-Jacobi equations and also discuss the computational efficiency of the algorithm. We then illustrate the inner-workings and nuances of our new stochastic DO level set energy optimal path planning algorithm through two simple, yet important, canonical steady flows that simulate a steady front and a steady eddy. We formulate a double energy-time minimization to obtain a semi-analytical energy optimal path for the steady front crossing test case and compare the results to these of our stochastic DO level set scheme. We then apply our methodology to an idealized ocean simulation using Double Gyre flows, and finally show an application with real ocean data for completing a mission in the Middle Atlantic Bight and New Jersey Shelf/Hudson Canyon region.

The ocean is a complex, constantly changing, highly dynamical system. Prediction capabilities are constantly being improved in order to better understand and forecast ocean properties for applications in science, industry, and maritime interests. Our overarching goal is to better predict the ocean environment in regions of complex topography with a continental shelf, shelfbreak, canyons and steep slopes using the MIT Multidisciplinary Simulation, Estimation and Assimilation Systems (MSEAS) primitive-equation ocean model. We did this by focusing on the complex region surrounding Taiwan, and the period of time immediately following the passage of Typhoon Morakot. This area and period were studied extensively as part of the intense observation period during August – September 2009 of the joint U.S. – Taiwan program Quantifying, Predicting, and Exploiting Uncertainty Department Research Initiative (QPE DRI). Typhoon Morakot brought an unprecedented amount of rainfall within a very short time period and in this research, we model and study the effects of this rainfall on Taiwan’s coastal oceans as a result of river discharge. We do this through the use of a river discharge model and a bulk river-ocean mixing model. We complete a sensitivity study of the primitive-equation ocean model simulations to the different parameters of these models. By varying the shape, size, and depth of the bulk mixing model footprint, and examining the resulting impacts on ocean salinity forecasts, we are able to determine an optimal combination of salinity relaxation factors for highest accuracy.

p> The science of uncertainty quantification has gained a lot of attention over recent years. This is because models of real processes always contain some elements of uncertainty, and also because real systems can be better described using stochastic components. Stochastic models can therefore be utilized to provide a most informative prediction of possible future states of the system. In light of the multiple scales, nonlinearities and uncertainties in ocean dynamics, stochastic models can be most useful to describe ocean systems.

Uncertainty quantification schemes developed in recent years include order reduction methods (e.g. proper orthogonal decomposition (POD)), error subspace statistical estimation (ESSE), polynomial chaos (PC) schemes and dynamically orthogonal (DO) field equations. In this thesis, we focus our attention on DO and various PC schemes for quantifying and predicting uncertainty in systems with external stochastic forcing. We develop and implement these schemes in a generic stochastic solver for a class of non-autonomous linear and nonlinear dynamical systems. This class of systems encapsulates most systems encountered in classic nonlinear dynamics and ocean modeling, including flows modeled by Navier-Stokes equations. We first study systems with uncertainty in input parameters (e.g. stochastic decay models and Kraichnan-Orszag system) and then with external stochastic forcing (autonomous and non-autonomous self-engineered nonlinear systems). For time-integration of system dynamics, stochastic numerical schemes of varied order are employed and compared. Using our generic stochastic solver, the Monte Carlo, DO and polynomial chaos schemes are intercompared in terms of accuracy of solution and computational cost.

To allow accurate time-integration of uncertainty due to external stochastic forcing, we also derive two novel PC schemes, namely, the reduced space KLgPC scheme and the modified TDgPC (MTDgPC) scheme. We utilize a set of numerical examples to show that the two new PC schemes and the DO scheme can integrate both additive and multiplicative stochastic forcing over significant time intervals. For the final example, we consider shallow water ocean surface waves and the modeling of these waves by deterministic dynamics and stochastic forcing components. Specifically, we time-integrate the Korteweg-de Vries (KdV) equation with external stochastic forcing, comparing the performance of the DO and Monte Carlo schemes. We find that the DO scheme is computationally efficient to integrate uncertainty in such systems with external stochastic forcing.

A new methodology for Bayesian inference of stochastic dynamical models is developed. The methodology leverages the dynamically orthogonal (DO) evolution equations for reduced-dimension uncertainty evolution and the Gaussian mixture model DO filtering algorithm for nonlinear reduced-dimension state variable inference to perform parallelized computation of marginal likelihoods for multiple candidate models, enabling efficient Bayesian update of model distributions. The methodology also employs reduced-dimension state augmentation to accommodate models featuring uncertain parameters. The methodology is applied successfully to two high-dimensional, nonlinear simulated fluid and ocean systems. Successful joint inference of an uncertain spatial geometry, one uncertain model parameter, and 0(105) uncertain state variables is achieved for the first. Successful joint inference of an uncertain stochastic dynamical equation and 0(105) uncertain state variables is achieved for the second. Extensions to adaptive modeling and adaptive sampling are discussed.

Data assimilation, as presented in this thesis, is the statistical merging of sparse observational data with computational models so as to optimally improve the probabilistic description of the field of interest, thereby reducing uncertainties. The centerpiece of this thesis is the introduction of a novel such scheme that overcomes prior shortcomings observed within the community. Adopting techniques prevalent in Machine Learning and Pattern Recognition, and building on the foundations of classical assimilation schemes, we introduce the GMM-DO filter: Data Assimilation with Gaussian mixture models using the Dynamically Orthogonal field equations.

We combine the use of Gaussian mixture models, the EM algorithm and the Bayesian Information Criterion to accurately approximate distributions based on Monte Carlo data in a framework that allows for efficient Bayesian inference. We give detailed descriptions of each of these techniques, supporting their application by recent literature. One novelty of the GMM-DO filter lies in coupling these concepts with an efficient representation of the evolving probabilistic description of the uncertain dynamical field: the Dynamically Orthogonal field equations. By limiting our attention to a dominant evolving stochastic subspace of the total state space, we bridge an important gap previously identified in the literature caused by the dimensionality of the state space.

We successfully apply the GMM-DO filter to two test cases: (1) the Double Well Diffusion Experiment and (2) the Sudden Expansion fluid flow. With the former, we prove the validity of utilizing Gaussian mixture models, the EM algorithm and the Bayesian Information Criterion in a dynamical systems setting. With the application of the GMM-DO filter to the two-dimensional Sudden Expansion fluid flow, we further show its applicability to realistic test cases of non-trivial dimensionality. The GMM-DO filter is shown to consistently capture and retain the far-from-Gaussian statistics that arise, both prior and posterior to the assimilation of data, resulting in its superior performance over contemporary filters. We present the GMM-DO filter as an efficient, data-driven assimilation scheme, focused on a dominant evolving stochastic subspace of the total state space, that respects nonlinear dynamics and captures non-Gaussian statistics, obviating the use of heuristic arguments.

A new generation of efficient parallel, multi-scale, and interdisciplinary ocean models is required for better understanding and accurate predictions. The purpose of this thesis is to quantitatively identify promising numerical methods that are suitable to such predictions. In order to fulfill this purpose, current efforts towards creating new ocean models are reviewed, an understanding of the most promising methods used by other researchers is developed, the most promising existing methods are studied and applied to idealized cases, new methods are incubated and evaluated by solving test problems, and important numerical issues related to efficiency are examined. The results of other research groups towards developing the second generation of ocean models are first reviewed. Next, the Discontinuous Galerkin (DG) method for solving advection-diffusion problems is described, including a discussion on schemes for solving higher order derivatives. The discrete formulation for advection-diffusion problems is detailed and implementation issues are discussed. The Hybrid Discon- tinuous Galerkin (HDG) Finite Element Method (FEM) is identified as a promising new numerical scheme for ocean simulations. For the first time, a DG FEM scheme is used to solve ocean biogeochemical advection-diffusion-reaction equations on a two- dimensional idealized domain, and p-adaptivity across constituents is examined. Each aspect of the numerical solution is examined separately, and p-adaptive strategies are explored. Finally, numerous solver-preconditioner combinations are benchmarked to identify an efficient solution method for inverting matrices, which is necessary for implicit time integration schemes. From our quantitative incubation of numerical schemes, a number of recommendations on the tools necessary to solve dynamical equations for multiscale ocean predictions are provided.

A fundamental requirement in realistic computational geophysical fluid dynamics is the optimal estimation of gridded fields and of spatial-temporal scales directly from the spatially irregular and multivariate data sets that are collected by varied instruments and sampling schemes. In this work, we derive and utilize new schemes for the mapping and dynamical inference of ocean fields in complex multiply-connected domains, study the computational properties of our new mapping schemes, and derive and investigate new schemes for adaptive estimation of spatial and temporal scales.

Objective Analysis (OA) is the statistical estimation of fields using the Bayesian- based Gauss-Markov theorem, i.e. the update step of the Kalman Filter. The existing multi-scale OA approach of the Multidisciplinary Simulation, Estimation and Assimilation System consists of the successive utilization of Kalman update steps, one for each scale and for each correlation across scales. In the present work, the approach is extended to field mapping in complex, multiply-connected, coastal regions and archipelagos. A reasonably accurate correlation function often requires an estimate of the distance between data and model points, without going across complex land- forms. New methods for OA based on estimating the length of optimal shortest sea paths using the Level Set Method (LSM) and Fast Marching Method (FMM) are derived, implemented and utilized in general idealized and realistic ocean cases. Our new methodologies could improve widely-used gridded databases such as the climatological gridded fields of the World Ocean Atlas (WOA) since these oceanic maps were computed without accounting for coastline constraints. A new FMM-based methodology for the estimation of absolute velocity under geostrophic balance in complicated domains is also outlined. Our new schemes are compared with other approaches, including the use of stochastically forced differential equations (SDE). We find that our FMM-based scheme for complex, multiply-connected, coastal regions is more efficient and accurate than the SDE approach. We also show that the field maps obtained using our FMM-based scheme do not require postprocessing (smoothing) of fields. The computational properties of the new mapping schemes are studied in detail. We find that higher-order schemes improve the accuracy of distance estimates. We also show that the covariance matrices we estimate are not necessarily positive definite because the Weiner Khinchin and Bochner relationships for positive definiteness are only valid for convex simply-connected domains. Several approaches to overcome this issue are discussed and qualitatively evaluated. The solutions we propose include introducing a small process noise or reducing the covariance matrix based on the dominant singular value decomposition. We have also developed and utilized novel methodologies for the adaptive estimation of spatial-temporal scales from irregularly spaced ocean data. The three novel methodologies are based on the use of structure functions, short term Fourier transform and second generation wavelets. To our knowledge, this is the first time that adaptive methodologies for the spatial-temporal scale estimation are proposed. The ultimate goal of all these methods would be to create maps of spatial and temporal scales that evolve as new ocean data are fed to the scheme. This would potentially be a significant advance to the ocean community for better understanding and sampling of ocean processes.

In this thesis, we conduct research toward understanding coupled physics-biology processes in ocean straits. Our focus is on new analytical studies and higher-order simulations of idealized dynamics that are relevant to generic biological processes. The details of coupled physics-biology models are reviewed and an in-depth global equilibrium and local stability analysis of a Nutrient-Phytoplankton-Zooplankton (NPZ) model is performed. This analysis includes parameter studies and methods to evaluate parameter sensitivity, especially in the case where some system parameters are unknown. As an initial step toward investigating the interaction between physics and biology in ocean straits, we develop and verify a new coupled physics-biology model for two-dimensional idealized physical processes including tides and apply it to the San Bernardino Strait in the Philippine Archipelago. This two-dimensional numerical model is created on a structured grid using operator splitting and masking. This model is able to accurately represent biology for various physical flows, including advection-dominated flows over discontinuities, by using the Weighted Essentially Non-Oscillatory (WENO) scheme. The numerical model is verified against a Discontinuous-Galerkin (DG) numerical scheme on an unstructured grid. Several simulations of tidal flow are completed using bathymetry and flow magnitudes com- parable to those found in the San Bernardino Strait with different sets of parameters, tidal periods, and levels of diffusion. Results are discussed and compared to those of a three-dimensional modeling system. New results include: new methods for analyzing stability, the robust two-dimensional model designed to best represent advection-dominant flows with minimal numerical diffusion and computational time, and a novel technique to initialize three-dimensional biology fields using satellite data. Additionally, application of the two-dimensional model with tidal forcing to the San Bernardino Strait reveals that flow frequencies have strong influence on biology, as very fast oscillations act to stabilize biology in the water column, while slower frequencies provide sufficient transport for increased biological activity.

In this thesis, we explore the different methods for parameter estimation in straightforward diffusion problems and develop ideas and distributed computational schemes for the automated evaluation of physical and numerical parameters of ocean models. This is one step of “adaptive modeling”. Adaptive modeling consists of the automated adjustment of self-evaluating models in order to best represent an observed system. In the case of dynamic parameterizations, self-modifying schemes are used to learn the correct model for a particular regime as the physics change and evolve in time.

The parameter estimation methods are tested and evaluated on one-dimensional tracer diffusion problems. Existing state estimation methods and new filters, such as the unscented transform Kalman filter, are utilized in carrying out parameter estimation. These include the popular Extended Kalman Filter (EKF), the Ensemble Kalman Filter (EnKF) and other ensemble methods such as Error Subspace Statistical Estimation (ESSE) and Ensemble Adjustment Kalman Filter (EAKF), and the Unscented Kalman Filter (UKF). Among the aforementioned recursive state estimation methods, the so-called “adjoint method” is also applied to this simple study.

Finally, real data is examined for the applicability of such schemes in real-time fore- casting using the MIT Multidisciplinary Simulation, Estimation, and Assimilation System (MSEAS). The MSEAS model currently contains the free surface hydrostatic primitive equation model from the Harvard Ocean Prediction System (HOPS), a barotropic tidal prediction scheme, and an objective analysis scheme, among other models and developing routines. The experiment chosen for this study is one which involved the Monterey Bay region off the coast of California in 2006 (MB06). Accurate vertical mixing parameterizations are essential in this well known upwelling region of the Pacific. In this realistic case, parallel computing will be utilized by scripting code runs in C-shell. The performance of the simulations with different parameters is evaluated quantitatively using Pattern Correlation Coefficient, Root Mean Squared error, and bias error. Comparisons quantitatively determined the most adequate model setup.