My overarching research goal is to bring the advanced tools of mathematics to bare on relevant real world problems bringing new perspectives to some fields, improving general scientific understanding, and fostering new methods for prediction and analysis. In my most recent work, I have been exploring the interplay between mathematics and data primarily through the lens of Data Assimilation (DA) and Statistical or Machine Learning techniques. Working with others from across the globe we have: developed new techniques allowing the application of the Ensemble Kalman Filter (EnKF) to models which use non-conservative adaptive mesh solvers [23]; applied ensemble smoothing techniques with an age stratified SEIR model to study the spread of SARS-CoV-2 (the virus that causes COVID-19) obtaining continuous time series of the basic reproductive number R₀ which, can be used to evaluate successes and failures of various mitigation strategies [6]; extended the previously mentioned model and DA method to a multi-group population model to study disparities in the impact of the pandemic on different racial/ethnic groups [8]; and employed new metrics using boundary preserving diffeomorphisims to calibrate a sea ice model [27] (that explicitly models leads) using SAR sattelite imagery [11]. I am also currently involved in projects aimed at: using a reduced order model of the polar vortex [19] in combination with reanalysis data from the European Center for Medium Weather Forecasting (ECMWF) to understand bi-stability in the vortex itself; employing Generative Adversarial Networks (GANs) trained on images of melt ponds on Arctic sea ice to produce realistic meltpond geometries for use in sea ice modeling; investigating how machine learning techniques can be used as an augmentative tool in data analysis and Data Assimilation; and studying sea ice microstructure through the lens of topological data analysis [21].

In my graduate studies, I have focused on the physical properties and processes of sea ice and their role in the Earth’s climate system. I have been involved in research which applies: percolation theory to electrical and fluid transport through the porous microstructure of sea ice [10],[9]; inversion theory to problems in remote sensing [20]; graph theory to melt pond evolution [1]; and homogenization theory to ocean wave dynamics in the marginal ice zones (MIZ) of the Arctic and Antarctic [22]. That work was been guided by the fact that sea ice may be thought of as a composite material over multiple scales. From the sub millimeter scale of its porous structure to the kilometer scale of the mix of ice floes and open water, the tools of homogenization theory consistently apply.

I have also been fortunate enough to complement my mathematical research with experimental field work in both the Arctic, Antarctic and Gulf of Alaska. This has given me a unique perspective on the challenges associated with combining theoretical work and experimental science. I have a deep interest in mathematics and its application and seek to continue my work bringing the two together.

In general, my mathematical interests include: data assimilation, data analysis, machine learning, fluid dynamics , continuum mechanics , dynamical systems, homogenization theory, inversion theory, numerical analysis, PDE’s, topological data analysis and complex analysis. My scientific interests include: climate science, polar and marine science, mathematical epidemiology, mathematics and social justice, statistical physics, and planetary science.

Data Assimilation and Ensemble Methods

Ensemble Methods for Parameter Estimation and General Reanalysis

Ensemble smoother techniques can be derived by assuming a perfect forward model.

(1)

In general, x is the realization of model parameters, and y consists of the uniquely predicted measurements resulting from the parameters x through the model operator g(x). The vector y then consists of the predicted measurements given some model error, e. To solve the inverse problem, it is efficient to frame it into an equation using Bayes’ theorem:

(2)

Equation (2) represents the so-called smoothing problem, which can be approximated using ensemble methods. In practice this problem can be solved using an iterative ensemble smoothing method called an ensemble smoothing with multiple data assimilation (ESMDA). The method, which is similar to an ensemble smoother, solves the parameter-estimation problem and is formally derived from the Bayesian formulation using a tempering procedure [26]. What separates this method from other similar methods is that it approximates the posterior recursively, gradually introducing information to alleviate the impact of any nonlinear approximations. After updating through to the last time step, it begins the assimilation process over again - resampling the vector of perturbed observations, which reduces sampling error [5]. The general steps can be described simply. (1) Sample a large ensemble of the parameters you wish to estimate assigning uncertainty to each. (2) Integrate the ensemble of model realizations forward in time to obtain a prior ensemble prediction. (3) Compute the posterior of the ensemble of parameters using the difference between predicted states and observations with the correlations between those parameters and predicted states. (4) Integrate the ensemble of updated model realizations forward in time and repeat for the number of ESMDA steps chosen. After the cycle is completed the ensemble mean of parameters are considered optimal as well as the mean of model predictions arising from the ensemble of parameters. One also obtains a measure of uncertainty through the spread in the parameters as well.

These methods have broad applicability and are only constrained by dimensionallity and the computational cost of the forward model. In addition to the advantage of providing a measure of uncertainty ESMDA can be used to estimate critically important parameters through the combination of the model and available data in the reanalysis sense. A particularly relevant example right now is that of current COVID-19 pandemic. In [7] a Susceptible Infected Exposed and Recovered (SEIR) model which was compartmentalized through an age stratified contact rate matrix used ESMDA to estimate important parameters and make short term predictions by a large group of international researchers in their regions. Participating in this effort our group used this model to compare the success of mitigation efforts of four U.S. states by estimating the basic reproduction number of the SARS-CoV-2 virus as a function of time R(t). Shown in Figure 1 is a sample of the results where a gradual step down of R(t) is used starting at dates which correspond to sufficient decreases in mobility data as measured from cell phones. In this figure we highlight that New York achieved a consistent R(t) < 1 implying decreasing cases while California fluctuated around R(t) = 1 during the data window of early March to May 20^th. The model was conditioned on cumulative deaths and daily hospitalizations. Methods like these can be used to preform a reanalysis of the spread of SARS-CoV-2 and tied to other historical data to determine successful strategies and better prepare for the future.

California	New York

Figure 1: US Case 3: A gradual step down in R(t) with the first and intermediate values chosen so that the prior mean closely follows the data until the time for which R(t) is guessed to be one. The red thin line in the plots for R(t) is an indication of the value R(t) = 1 for easier identification.

The model in [7] has also recently been extended to model sub-populations with transmission between populations possible as well as the provision of a continuous prior on R(t). The sub-populations can represent different countries, states, other localities or population groups with in a locality. In [8] we use cumulative death data reported for different racial/ethnic groups to analyze and study observed disparities in the spread of SARS-CoV-2 across these groups in the United States. We focus on 9 states with the most complete reporting of this data as well as the District of Columbia.

Case: Continuous Updates to R_t(n) (HI)

Figure 2: Analysis results for the continuous update case for the state of District of Columbia .

In this more recent analysis a continuous prior for the given locality is taken from rt.live and used for all groups and updated based on the death data over time. In Figure ?? we see the results for the District of Columbia which shows that after the mitigation period (corresponding to the large dip in R(t)) a stratification in the basic reproduction number occurs across the different racial/ethnic groups with the Black and Latinx communities showing the highest rate of spread and more deaths than other groups. This is particularly notable since the Black community is roughly the same proportion of the population of D.C. while the Latinx community only about 11%. We see this as the detection of front-line communities who bores the brunt of the pandemic for a variety of reasons. Methods like ESMDA can be used in this way to shed light on inequity, determine disparity and provide social scientists and policy makers with valuable estimates of critical parameters needed for decision making and equitable actions.

Future Directions: These methods are widely applicable to a broad range of problems, currently I am involved in using methods like this to investigate bi-stability in the polar vortex through the lens of a reduced order model and data from The ECMWF. I also envision extending the kind of analysis discussed above to continue studies of the COVID-19 pandemic in the context of social justice and assessment of various mitigation measures. The ESMDA method is model agnostic at its heart and more complex models could also be employed i.e. agent based or network models. This type of data assimilation is also useful for short term forecasting as any number of model parameters can be estimated including case fatality rates, reproduction numbers etc. As a result a system which estimates parameters that can serve as early warning signs of trouble could be used during ongoing outbreaks in the future for rapid response. In addition, I believe it possible to develop a general model of hospitalization rates, for which a general epidemiological model of viral spread is a part, and through the assimilation of those rates the early detection of a new outbreak. I am also interested in studying how climate change had disproportionately affected marginalized groups and see methods like this also useful. Air quality modeling in conjunction with air quality and health data, storm surge and flood models, extreme temperature events and vulnerable populations in conjunction with historical data are just some of the other topics that can be explored. There is also the potential of methods like these to provide undergraduate students with valuable research experience as general codes can be given to them allowing for accessible exploration of a specific model, data assimilation methodology, and relevant data set. The broad application of these methods also provides for broad appeal and can facilitate mentoring oppurtunities for students from a variety of backgrounds.

Ensemble Kalman Filter for Adaptive Moving Mesh Models

Ensemble data assimilation has been widely implemented in weather prediction [12], tumor growth and spread [15], and petroleum reservoir history matching [25]. The ensemble Kalman filter (EnKF) relies on estimates of error statistics using an ensemble of model runs assumed to be Gaussian distributed. The error estimates themselves are calculated using the state vectors (x_i^f) formed from the state variables of each ensemble member. When observations are available each of the N^e ensemble members are updated according to

(3)

Here, y is a vector of observations, h is the operation operator which maps model variables to the observations and K is the Kalman gain matrix and has the form,

(4)

The i^th column of the matrix X^f is formed from the difference between each state vector x_i^f and their mean x^f and normalized by the number of ensemble members minus one. The i^th column matrix Y is formed from the difference between each state vector and mean after applying the observation operator h to each and R^e is the observation error covariance matrix and can also be estimated from the ensemble statistics. The EnKF is often used in a sequential mode by applying updates at times when observations are available and using model to forecast in between.

The EnKF can be directly applied to any forward in time model for which all ensemble members have the same dimension however some adaptations may be needed in special cases such as models which use adaptive mesh solvers. Numerical solvers using adaptive meshes can focus computational power on important regions of a model domain capturing important or unresolved physics. The adaptation can be informed by the model state, external information, or made to depend on the model physics. This has caused AMM schemes to rise in popularity and thus adaptation of DA schemes like the EnKF become important and even advantageous. In recent work we developed an EnKF scheme for a 1-d adaptive mesh solver that in which the nodes are advected with the flow and sometimes deleted if two nodes are too close together (with in a tolerance of δ₁) or inserted if two nodes are two far apart (outside a tolerance of δ₂). While AMM schemes like this are useful to focus computational power in regions of large gradients they present significant challenges for ensemble DA methods. This is because each ensemble member may have a different number of nodes all in different locations making the calculation of the needed error statistics difficult. The different number of nodes can be addressed by interpolating to new nodes for each ensemble member to match the dimension across the ensemble. However, this does not account for the fact that nodes are in different locations and more care is needed. There are two primary options for addressing this, (1) interpolating each ensemble member to the same grid or (2) including the node locations in the state vector paring nodes in some reasonable way. For the scheme described above we use a partitioning of the domain based on the closeness tolerance δ₁. Partitioning the domain into intervals of size δ₁ there is a guarantee that each interval will have at most one point with in it. From there one can map each ensemble member to reference mesh defined by the partition or interpolate new nodes for each ensemble member that does not have a point in a given interval. In the former case the update is done on a regular grid with all nodes in the same location and each ensemble member can be mapped back to the mesh structure it had before the update. In the latter case we include the node locations in the state vector comparing nodes which occupy the same intervals to calculate the needed error statistics for the EnKF. In this way the node locations themselves are updated based on any available physical observations, which are observations of the quantities that drive that flow. While the larger state vector does increase computational cost we did find that the inclusion and update of those node locations improved the assimilation across several metrics including the standard RMSE measure compared to instead only updating the physical quantities on the reference mesh defined by the tolerance δ₁. The models we consider in that work are Burgers equation and the Kuramoto Sivashinsky equation. In Figure 3 we show the blocks of the error covariance matrices, calculated as X^f(X^f)^T, comparing covariance between the physical values and the node locations themselves. Strikingly we see that the diagonal of that block takes on a very similar shape as that of the first spatial gradient of the solutions highlighting the extra information brought into the update.

Future Directions: While the work described was in 1-d, the complexity and range of patterns that can form in 2-d or 3-d is far greater than is possible in 1-d. This implies that significantly more information may be carried in the cross-covariances between physical values and node locations. Further, if the motivation for the use of an AMM scheme is to focus computational power in regions of strong gradients, updating those node locations in accordance with where observations of those gradients are large may be very advantageous. If the remeshing rules for the AMM model are based on strict considerations of node distances and mesh geometries, a 2-d or 3-d analogue of the reference mesh should be attainable enabling the application of a scheme which includes the mesh in the . An example of such a model is the novel Lagrangian sea ice model neXtSIM [18, 3] which uses a non-conservative adaptive triangular mesh and finite elements. The key to implementing a scheme which updates the mesh is the configuration of an augmented vector that includes variables characterizing the underlying numerical mesh. In 1-d there is little ambiguity as the natural approach is to just add the grid points of the underlying mesh. But in 2-d, the characterization of the “grid in terms of such vector is not as straightforward when, for instance, using a finite element method. In that case Physical values may be stored at the centers of the elements, where as the elements themselves are determined by their vertices. Considerations of the mesh type, triangular, hexagonal, square etc. also add complications. The development of 2-d or 3-d mesh update schemes is a rich mathematical problem involving numerical analysis, statistics and topology. It also promises to relate to a wide range of applications.

New Metrics for Model Calibration

For models which produce coherent structures such as ocean, atmospheric, or sea ice models, comparison of model output to observational data through something like a standard euclidean metric might be misleading or not informative enough. In the presence of coherent structures such as hurricanes, eddies, or sea ice leads one may want to measure model skill with considerations of the geometries as opposed to a simple point by point comparison of values at specific locations. Using the Elastic-Dechohesive rheology the model presented in [27] explicitly represents sea ice leads through discontinuities tracked by a jump vector. While predicting the location of leads is important the shape and orientation of the cracks is also an important consideration and warrants a metric which can take this into account. This can be accomplished by warping the model out put, viewed as an image, to observational data first optimally aligning the leads and then calculating a difference metric. The warping is done with a boundary preserving diffeomorphism and provides two metrics: how close two images can possibly be and how much warping it took to warp one image to another. The space of images considered is the set F = {f : D → Rⁿ|f ∈ C^∞(D)} with dim(D) = m ≤ n they are warped using the set of orientation preserving diffeomorphisims γ ∈ Γ that preserve the boundary of D which forms a group under the operation of composition f ∘ γ. The metric which gives the distance between two images f₁,f₂ ∈ F is defined as,

(5)

This distance, the amplitude distance, satisfies all of the requirements of a proper metric and measures how close two images are after they are optimally aligned by the diffeomorphism γ the amount of warp from this gamma, the phase distance, can also be used as a measure of how different the two images are. These metrics were applied to model out put for the model in [27] compared to SAR satellite imagery for a variety of shear and tensile strength parameters sampled using a latin hyper cube sampling scheme. Using the metrics above a statistical emulator was employed to estimate the optimal parameters based on the data. In Figure 4 we can see that when using the two metrics above a meaningful calibration was obtained while simply using the Euclidean distance provided a flat posterior with little new information. All of the details can be found in [11]

Figure 4: Posteriors for model calibration. Left: Using the amplitude and phase distances. Right: Using just the Euclidean distance.

Future Directions: Metrics like these can be applied to a variety of models which produce coherent structures directly. However, finding the optimal gamma is difficult and new optimization techniques, rather than just gradient dissent are desirable. Improving the optimization algorithm is a part of on going work. I also have deep interest in the development of new metrics and studies of which metrics fit a specific situation best.

Machine Learning, Sea Ice and Geophysics

Machine Learning (ML) has become a hot topic in applied mathematics and statistics in recent years owing to the direct utility some of the methods allow. I see ML as powerful for data analysis and an augmentative tool for modeling and forecast. I am currently involved in two projects taking advantage of some of promise of ML techniques to aid in sea ice modeling efforts and as an augmentative tool for Data Assimilation.

In the first we are exploring the use of a Generative Adversarial Network (GAN) trained on images of melt ponds atop sea ice to generate realistic melt pond geometries and scenes which can serve as input to large scale sea ice models for which parameters may depend on sub-scale processes related to the ponds themselves. For example, ice strength is directly tied to ice thickness in most numerical sea ice models and melt ponds have an effect on the ice thickness over the area which they cover. If, based on thermodynamic variables in the model, the presence of ponds is expected one may wish to generate a statistically consistent melt pond geometry to inform the strength parameter in the ice model. This is particularly important for models wishing to use large ensembles as many realistic melt pond scenes may be needed. Some preliminary results are shown in Figure  5.

GAN Generated	Real Binarized	Area v.s. Perimeter

Figure 5: Examples of GAN generated vs Real (but binarized) images of melt pond scenes with a comparision of the Area vs Perimeter relationship for the generated set (Blue) and Real (Red) ponds.

An additional use for such GAN could be in the aid of sea ice concentration retrievals. In the Arctic summer months melt ponds cover vast areas of sea ice which have the same microwave signature as that of open water. Passive microwave radiometry uses the contrast in microwave signatures of sea ice and open water to estimate sea ice concentration usually at or above the 14km scale. When the melt ponds sit atop the ice they obscure the ice beneath which can cause up to a 30% in sea ice concentration retrieval [14]. If one can rapidly generate sea ice scenes with open water and melt ponds through a well trained GAN, which, are consistent with specific satellite observed microwave radiance. Then it may possible to obtain a measure of uncertainty in what concentrations actually produced the signal. Further, using available validated data it may be possible to train a Deep Neural Net to select the most likely candidate of those possibilities.

The second project looks using ML to replace the observation operator in Data Assimilation. This is motivated by the sea ice concentration example but is a generally applicable idea. The observation operator maps model state variables to the observations space where the predictions and observations may be compared. For passive microwave retrieved sea ice concentrations one may choose to use the concentration value provided by NASA ( with up to 30% due to ponding), for example in which case the observation operator is a simple projection or one may map the sea ice state what would be the observable microwave radiances. The former case is computationally inexpensive and the latter not, especially when using ensemble methods. Assimilating data that may be in in 30% error is undesirable and can be avoided by instead mapping sea ice state to satellite radiance. ML enables this through cuts in computational cost. This concept could apply to many other geophysical systems of interest.

Future Directions: In addition to the projects described above I am also interested in finding ways to use ML to augment models and data assimilation schemes. For computationally expensive models the ability to run large ensembles is often a missing but critical factor in improving forecast skill. This is particularly true with ensemble DA methods where larger ensemble sizes provide better error statistics. While ML can not yet fully replace a numerical model in long term forecasting there has been much progress in providing some short term prediction skill [16, 17]. In the ensemble DA setting this is sufficient as more ensemble members can be created a short time before an observation is available and progressed to the observation time using the ML forecasting scheme. This could potentially allow the use of ensemble methods in cases where computational cost is prohibitive as such. I am also interested in using self learning techniques such a self organizing maps in data analysis in general. ML and its interplay with applied mathematics presents many opportunities.

Wave Dynamics in the Marginal Ice Zone

The MIZ can be loosely defined as the region of the ice pack with an ice concentration low enough such that it is significantly effected by ocean swell. Waves propagating into the ice pack play a central role in determining the the break up, distribution, and sizes of ice floes in the MIZ. The ice floes in turn control the dispersion and attenuation of the of the incoming waves. These waves are important to consider in any sea ice model and play a major role in ice melt, shore erosion in the Arctic, and the polar ecosystem. Recently, continuum models have been developed which treat the MIZ as a two-component composite of ice and slushy water [28, 13]. More specifically, the ice-water composite has been modeled as a thin elastic plate, a viscous ice layer atop an inviscid fluid, as well as a viscoelastic layer atop an in-viscid fluid. These models yield dispersion relations which determine the wave properties necessary for continued propagation through the MIZ. At the heart of these models are effective parameters, namely, the effective elasticity, viscosity, and complex viscosity. In practice, these effective parameters, which depend on the composite geometry and the physical properties of the constituents, are quite difficult to measure. To overcome this difficulty, we have employed the tools of homogenization theory to derive bounds on the effective parameters necessary for these models [22]. Specifically, we have developed a Stieltjes integral representation, involving a positive spectral measure of a self-adjoint operator, for the effective complex viscoelasticity ν^* of the ice water composite in the low frequency, long wave limit and have the form,

(6)

(7)

The positivity of the measure allows for the derivation of rigorous forward bounds which depend on partial information about the geometry of the ice-water composite. Integral representations like this have also been used to derive inverse bounds, where measured data is used to bound the structural parameters of the composite, such as the volume factions of the components [4]. We have also developed a simplified wave equation for waves in the ice-water composite. This wave equation arises naturally from the homogenization of the governing equations in the low frequency limit, and could provide a computationally inexpensive numerical model of wave propagation in the MIZ.

General Bounds	Matrix Particle Bounds

Figure 6: Left: General bounds derived from the Stieltjes integral representations for various wave periods T. Right: Tighter bounds using the restricted support under the matrix particle assumption for T = 25

The bounds are shown in Figure 6 as well as an extension of the theory to the matrix-particle bounds valid in pancake ice. Pancake ice consists of small rounded ice floes embedded in a slushy layer typically separated by some small amount. With this particular geometric assumption the support of the measure in the above integral representations may be restricted,

(8)

allowing for a tightening of the original bounds [24, 2]. Here the restriction (s_m,S_M) is governed by the separation of the particles embedded in the matrix. In this case for reasonably chosen individual viscoelasticites of each separate component we obtain much tighter bounds which fall on a data point measured from lab grown pancake ice in a wave tank [29].

Future Directions: Once wave propagation characteristics are understood for a given configuration of ice, the next question is, how does the propagation of the wave affect the configuration? Waves of sufficient amplitude can certainly break up larger ice floes, this changes the geometry and thus the wave propagation characteristics in turn. I am interested in studying this problem through the development of a numerical model of wave propagation and ice break up using some of the results outlined above. The effect that wave propagation has on the floe size and thickness distribution of sea ice is important to consider in any large scale sea ice model. When a floe is broken up by wave action it becomes more susceptible to lateral melting, advection and rafting.

References