Partial differential equations (PDE) describe a wide range of phenomena across science and engineering. Of particular interest in this project is a suite of geophysical problems, including (1) rainwater runoff over land surfaces, (2) meltwater flow through snowpack, and (3) the flow of glaciers. The complexity of these problems is far beyond the ability of classical techniques to give exact solutions, and computer simulations are required to understand the models. However, existing numerical methods fail to capture certain essential features present in the physical and mathematical systems. Frequently, these features take the form of inequality constraints -- glacier thickness or water depths cannot be negative. This project will develop new methods that provide accurate approximations to such models while fully respecting physical bounds constraints. Through the co-investigator's involvement in the Juneau Icefield Research Program, findings of this project will directly inform field research. The mathematical and computational techniques developed in this project will be incorporated in open source software projects that are widely used in research and educational endeavors, pushing forward the state of the art in scientific simulation. This project will train a doctoral student in mathematics and provide for undergraduate research experiences through the McNair Scholars program at Baylor University.<br/> <br/>Many common variational problems obey a maximum principle, but only a restricted class of numerical methods preserves this important feature. An alternative is to explicitly enforce the maximum principle by recasting the problem as a variational inequality subject to physical bounds constraints. Additionally, many important problems inherently arise as variational inequalities. This project develops techniques for bounds-constrained solution of partial differential equations and variational inequalities. These techniques will be developed with suite of challenging geophysical problems in mind. These target applications include meltwater flow through snow, rainwater runoff over land surfaces, and glacier dynamics. These models are nonlinear, time-dependent, coupled partial differential equations with positivity constraints on a thickness variable. At the discrete level, bounds constraints are applied on the Bernstein control net of the finite element space. Since this allows for high-order spatial accuracy, it is also important to obtain high temporal accuracy for dynamic problems. This project will adapt Runge-Kutta methods to ensure bounds constraints hold over time. For single-stage methods, one can advance the solution by posing a variational inequality that minimizes a defect subject to bounds constraints. For higher-order time-stepping methods based on polynomial approximation, such as Galerkin-in-time or collocation-type Runge-Kutta methods, one can force the collocating polynomial to satisfy the bounds constraints uniformly in time, which may be especially important for tightly coupled processes. In addition to method development and computational applications, research will address how to adapt stability and convergence theory for Runge-Kutta methods to the bounds-constrained setting. To facilitate broader adoption by the computational mathematics and geophysics communities, newly-developed methodology will be included in the Firedrake project and the icepack library for glacier modeling built on top of it. Aspects of this research will be included in applied mathematics graduate classes at Baylor University and in earth and space sciences at the University of Washington.<br/><br/>This award reflects NSF's statutory mission and has been deemed worthy of support through evaluation using the Foundation's intellectual merit and broader impacts review criteria.