Efficient optimization of a regional water elevation model with an automatically generated adjoint

Abstract

Calibration of unknown model parameters is a common task in many ocean model applications. We present an adjoint-based optimization of an unstructured mesh shallow water model for the Baltic Sea. Spatially varying bottom friction parameter is tuned to minimize the misfit with respect to tide gauge sea surface height (SSH) observations. A key benefit of adjoint-based optimization is that computational cost does not depend on the number of unknown variables. Adjoint models are, however, typically very laborious to implement. In this work, we leverage a domain specific language framework in which the discrete adjoint model can be obtained automatically. The adjoint model is both exactly compatible with the discrete forward model and computationally efficient. A gradient-based quasi-Newton method is used to minimize the misfit. Optimizing spatially-variable parameters is typically an under-determined problem and can lead to over-fitting. We employ Hessian-based regularization to penalize the spatial curvature of the friction field to overcome this problem. The SSH dynamics in the Baltic Sea are simulated for a 3-month period. Optimization of the bottom friction parameter results in significant improvement of the model performance. The results are especially encouraging in the complex Danish Straits region, highlighting the benefit of unstructured meshes. Domain specific language frameworks enable automated model analysis and provide easy access to adjoint modeling. Our application shows that this capability can be enabled with few efforts, and the optimization procedure is robust and computationally efficient.

\draftfalse\journalname

Journal of Advances in Modeling Earth Systems (JAMES)

Finnish Meteorological Institute, Helsinki, Finland Imperial College, London, United Kingdom

\correspondingauthor

Tuomas Kärnä[email protected]

{keypoints}

Adjoint-based optimization is used to optimize bottom friction coefficient in 2D water elevation model for the Baltic Sea

The discrete adjoint model is automatically generated by leveraging a symbolic representation of the discrete forward model equations

The optimization method is robust and results in significant improvement in the sea surface height performance at tide gauge locations

Plain Language Summary

Ocean circulation models have several unknown parameters that must be tuned for each application in order to produce physically meaningful results. The tuning process can be a very laborious and time consuming task. In this paper, we investigate an automated way to tune the model’s friction at the sea bed to minimize the model’s error in predicted sea surface height. The method is based on a novel way of defining the model’s equations which enables solving such optimization problems automatically. The methodology is tested in the Baltic Sea. The modeled sea surface height is compared against observations at several tide gauges. We show that by optimizing the bottom friction, the model’s capability to predict sea surface height improves significantly. Moreover, we show that the optimization process is robust and computationally efficient.

1 Introduction

Numerical modeling is an indispensable tool in oceanography and climate sciences. One of its major bottlenecks, however, is fitting the model state to observations. Model configuration must be carefully tuned and calibrated to replicate observational data. In addition to merely running simulations, one typically also needs to quantify the uncertainty of model predictions [Kalmikov \BBA Heimbach (\APACyear2014), Loose \BBA Heimbach (\APACyear2021)], or estimate unknown parameters [Heemink \BOthers. (\APACyear2002), Zaron \BOthers. (\APACyear2011), Zhang \BOthers. (\APACyear2011), Almeida \BOthers. (\APACyear2018), Warder \BOthers. (\APACyear2022)]. The bottom friction coefficient, for example, is commonly regarded as an unknown free parameter in coastal applications (e.g., \citeAzhang2011). In addition, as the amount of high-resolution observational data increases, sophisticated data assimilation methods are needed to synthesize the observation data and fill in the gaps with a physically meaningful ocean state. All of these use cases are examples of inverse modeling.

Adjoint models provide an efficient way to solve inverse problems. In short, adjoint models allow evaluating the gradient of model outputs with respect to the model’s internal state, forcing fields, or parameters. The drawback is that implementing the adjoint model is technically challenging and labor intensive [Marotzke \BOthers. (\APACyear1999), Heemink \BOthers. (\APACyear2002), Vidard \BOthers. (\APACyear2015)]. The adjoints of the ROMS and NEMO ocean models, for example, have been implemented manually, by differentiating each operator of the model [Moore \BOthers. (\APACyear2004), Vidard \BOthers. (\APACyear2015)]. Maintaining such a hard-coded adjoint, however, requires constant human intervention as the models evolve; both ROMS and NEMO adjoint models now lag behind the latest forward model versions. In the case of the MITgcm model, the adjoint has been derived via Automatic Differentiation (AD), i.e. by automatically differentiating the Fortran source code [Marotzke \BOthers. (\APACyear1999)]. The MITgcm adjoint has been used in a wide variety of applications, but despite its success, the AD-generated adjoint approach has not been widely adopted in other models. Therefore, for the majority of ocean models, adjoint capability is not available and inverse problems are difficult to solve, e.g. one has to rely on ensemble runs to obtain statistical estimates of the model’s sensitivity. In this paper, we demonstrate that by leveraging suitable high-level abstractions and code generation, the adjoint model can be generated and solved automatically, with very few efforts once the forward model has been implemented.

Adjoint models come in two flavors, continuous and discrete adjoints. In the former, the adjoint equations are derived directly from the continuous model equations and then discretized with a method of choice (e.g. finite volume (FV) or finite element (FE) method). In the discrete adjoint case, on the other hand, one differentiates the discrete (e.g. FV or FE) equations to derive an adjoint that is exactly compatible with the discrete forward model. Differentiating the model source code, as in the case of AD, also results in a discrete adjoint. There are advantages and disadvantages to both continuous and discrete adjoint approaches [Sirkes \BBA Tziperman (\APACyear1997)]. For example, the continuous adjoint approach is flexible in the sense that the user is able to choose different discretization methods for the forward and adjoint equations, but comes with the burden of being tedious and error-prone to implement. The discrete adjoint approach, on the other hand, is less flexible, but requires less user effort and produces gradients that are exactly compatible with the discretized forward model. This is greatly beneficial in applications – such as the one considered in this paper – where gradient-based optimization methods are used, since the convergence of such methods is ensured.

In this work, we use an unstructured mesh FE ocean model, Thetis [Kärnä \BOthers. (\APACyear2018)], for which the discrete adjoint model can be derived automatically. Thetis has been implemented in the generic Firedrake FE modeling framework [Rathgeber \BOthers. (\APACyear2016)]. Firedrake uses a domain specific language (Unified Form Language, [Alnæs \BOthers. (\APACyear2014)]) to describe the FE weak forms. An automated code generator [Homolya \BOthers. (\APACyear2018)] is then used to generate C code to evaluate the terms of the weak form. The equations are assembled as a linear system and solved in parallel with the PETSc solver library [Balay \BOthers. (\APACyear2021)]. Leveraging the symbolic representation, it is possible to derive the discrete adjoint model automatically by differentiating the symbolic FE equations [Farrell \BOthers. (\APACyear2013)]. The pyadjoint library handles the taping of the operator calls, solving the adjoint equation and evaluating the gradient of the cost function. Finally, a gradient-based quasi-Newton method is used to solve the optimization problem.

We use a discontinuous Galerkin FE 2D shallow water model, implemented within the Thetis framework, to simulate tidal and atmospherically-driven water elevation dynamics in the Baltic Sea. Our primary focus is on the Baltic Sea and the Danish Straits, but as the North Sea is tightly coupled to the dynamics, the two seas must be simulated as a single dynamical system [Daewel \BBA Schrum (\APACyear2013), Hordoir \BOthers. (\APACyear2019), Kärnä \BOthers. (\APACyear2021)]. We use the adjoint model to optimize the model’s bottom friction coefficient to improve the representation of water elevation dynamics with respect to tide gauge observations.

The two seas exhibit quite different water elevation dynamics. The North Sea is a tidal system with dominant semi-diurnal tides [Huthnance (\APACyear1991)]. Tidal range varies from over 6 m in the English channel to roughly 0.4 m in Skagerrak. Tides are effectively filtered out in the Danish waters and only weak tides ( $<$ 10 cm range) are observed in the Baltic Sea. The atmospherically-induced water elevation gradient across the Danish Straits controls the water volume exchange to and from the Baltic Sea in synoptic and longer time scales [Mohrholz \BOthers. (\APACyear2015), Gräwe \BOthers. (\APACyear2015), Mohrholz (\APACyear2018)]. Consequently, water elevation in the Baltic Sea is mainly governed by episodic filling and emptying of the basin, storm surges, as well as atmospherically-generated seiche oscillations whose typical periods range from 17 to 31 h [Leppäranta \BBA Myrberg (\APACyear2009)].

Previously, adjoint models have been used in several coastal ocean applications. \citeAmassmann2010 used an AD-based adjoint of an unstructured mesh finite element model to study the sensitivity of a North Sea water elevation model versus the bottom friction coefficient, bathymetry, and open boundary forcing. They found that the model skill was most sensitive to bottom friction, although in certain regions there was some overlap with sensitivity to bathymetry. Similarly, \citeAwarder2021 studied the spatial and temporal sensitivies of a storm surge model of the North Sea with respect to bottom friction, bathymetry and wind forcing using Thetis.

\citeA

almeida2018 use the Delft3D-FLOW model to simulate a section of the Columbia River. A continous adjoint was implemented for the model and used to optimize the model’s bathymetry based on surface velocity measurements. As the continuous adjoint model does not match exactly with the discrete forward model, a careful solution strategy was devised to avoid numerical instabilities.

\citeA

heemink2002 presented an adjoint-based inverse model of a 3D hydrostatic model and used it to optimize bottom friction, vertical viscosity coefficient and bathymetry in an application to the European Continental Shelf. The adjoint model was implemented manually. A common problem in parameter estimation is that the problem is typically underdetermined: in the case of a spatially-varying parameter field, there are far too many degrees of freedom compared to the amount of observation data. \citeAheemink2002 reduced the dimensionality of the problem by allowing the parameter fields to vary only in predefined subdomains of the model grid. The choice of subdomains was informed by inspecting the gradient of the cost function, calculated with the adjoint, and also expert knowledge of the system dynamics (e.g. amphidromic points).

\citeA

zhang2011 used a shallow water adjoint model to optimize constant and spatially-variable bottom friction coefficient for the Bohai Sea and Yellow Sea. Also here, the dimensionality of the problem was reduced by optimizing the friction values only in a sparse subset of the model grid points and using linear interpolation in between.

The availability of adjoint-based gradient information has a large number of other potentially powerful applications. In coastal engineering applications the adjoint can be used to optimize the design of engineering structures and their interaction with the environment. As an example, the aforementioned pyadjoint approach has been used to study the optimal placement of a large number of tidal turbines in [Funke \BOthers. (\APACyear2014), Funke \BOthers. (\APACyear2016), Goss \BOthers. (\APACyear2021)]. The automated adjoint approach can easily be extended to coupled model applications – as demonstrated in [Clare \BOthers. (\APACyear2022)] where inversion and sensitivity analyses where performed with Thetis coupled to a morphodynamic model. Finally, \citeAwallwork2020 demonstrated that adjoint-based goal-oriented error estimate can be used for mesh adaptation to optimize the model accuracy for a specified quantity of interest (e.g. power output of tidal turbines).

Compared to previous adjoint-based friction optimization applications, the novelty of the present work is the use of the automatically derived adjoint. This approach offers significant benefits to the user, as one only needs to implement the forward model and the cost function. This capability is built-in in Firedrake and therefore immediately available to any Firedrake-based implementation. In addition, we have implemented additional functionality in Thetis to provide the observation-based cost function, adjoint model related file IO – and the final missing ingredient – coupling with a quasi-Newton optimization method.

We use regularization to constrain the under-determined bottom friction optimization problem. That is, an additional term is added in the cost function to penalize the second derivatives of the friction field to control its smoothness. We show that the optimization procedure works well in a large and complex coastal domain where SSH dynamics vary across several time scales and the adjoint solver is computationally efficient.

The 2D shallow water model and its discretization is presented in Section 2. The variational inverse problem and its application to bottom friction optimization is presented in Section 3. Section 4 outlines the model configuration for the North Sea-Baltic Sea, followed by results in Section 5. Discussion and conclusions are presented in Sections 6 and 7, respectively. Further details on the forward and inverse modelling techniques are provided in A and B.

2 Shallow-water model

Denoting the free surface elevation and depth averaged velocity by $\eta$ and $\textbf{u}=(u,v)$ , respectively, the shallow water equations in non-conservative form read

	$\displaystyle\frac{\partial\eta}{\partial t}+\nabla\cdot(H\textbf{u})$	$\displaystyle=$	$\displaystyle 0$		(1)
	$\displaystyle\frac{\partial\textbf{u}}{\partial t}+\textbf{u}\cdot\nabla\textbf{u}+f\textbf{e}_{z}\wedge\textbf{u}+g\nabla\eta+\frac{1}{\rho_{0}}\nabla p_{a}$	$\displaystyle=$	$\displaystyle\textbf{D}_{\textbf{u}}+\frac{\bm{\tau}_{w}+\bm{\tau}_{b}}{H\rho_{0}},$		(2)

where $H=\eta+h$ is the total water column depth, $h$ is the bathymetry, $f$ is the Coriolis parameter, $\textbf{e}_{z}$ an upward vertical unit vector, $g$ the gravitational acceleration, $\rho_{0}$ the constant reference water density, $p_{a}$ the atmospheric pressure, $\textbf{D}_{\textbf{u}}=\nabla\cdot\big{(}\nu(\nabla\textbf{u}+(\nabla\textbf{u})^{T})\big{)}$ is the viscosity operator and $\nu$ is the horizontal eddy viscosity, $\bm{\tau}_{w}$ and $\bm{\tau}_{b}$ denote the surface (wind) and bottom stresses, respectively.

In this work, we use the wind stress formula by \citeAlarge2008:

	$\displaystyle\bm{\tau}_{w}$	$\displaystyle=$	$\displaystyle C_{D}^{w}\|\textbf{u}_{w}\|\textbf{u}_{w},$		(3)
	$\displaystyle C_{D}^{w}$	$\displaystyle=$	$\displaystyle\left\{\begin{array}[]{rl}a_{1}\frac{1}{\|\textbf{u}_{w}\|}+a_{2}+a_{3}\|\textbf{u}_{w}\|+a_{8}\|\textbf{u}_{w}\|^{6},&\|\textbf{u}_{w}\|<33\ \mathrm{m/s},\\ 0.00234,&\mathrm{otherwise}\end{array}\right.$		(6)

where $\textbf{u}_{w}$ stands for the wind velocity at 10 m height; the $a_{i}$ coefficients are defined in \citeAlarge2008.

The bottom friction is parametrized by the Manning formula,

	$\displaystyle\bm{\tau}_{b}$	$\displaystyle=$	$\displaystyle C_{D}\|\textbf{u}\|\textbf{u},$		(7)
	$\displaystyle C_{D}$	$\displaystyle=$	$\displaystyle g\frac{\mu^{2}}{H^{1/3}},$		(8)

where $\mu$ is the Manning friction coefficient. Generally $\mu$ can be regarded as an unknown, spatially-variable, model and mesh-dependent parameter.

2.1 Finite element discretization

We use the Thetis model implementation of the shallow water equations [Kärnä \BOthers. (\APACyear2018)]. The 2D model domain and its boundary are denoted by $\Omega$ and $\Gamma$ , respectively. The boundary consists of open ocean ( $\Gamma_{o}$ ) and closed land ( $\Gamma_{c}$ ) parts, i.e. $\Gamma=\Gamma_{o}\cup\Gamma_{c}$ .

Refer to caption — Figure 1: Finite elements: (a) linear continuous; (b) linear discontinuous; and (c) quadratic continuous scalar element. The dots denote scalar degrees of freedom. For vector-valued fields, each component is stored as a scalar.

Thetis supports several finite element families. In this work, the equations are discretized with linear discontinuous (DG₁) finite elements (see Fig. 1). The same element is used for both of the prognostic fields, $\eta$ and u. Let $\phi$ and $\bm{\phi}$ denote the scalar and vector valued test functions in the DG₁ function space. The weak form of the shallow water problem is obtained by multiplying the governing equations (1-2) by the test functions $\phi$ and $\bm{\phi}$ , respectively and integrating over the domain $\Omega$ :

	$\displaystyle\int_{\Omega}\frac{\partial\eta}{\partial t}\phi\;\mathrm{d}x+\mathcal{S}_{div}$	$\displaystyle=$	$\displaystyle 0,\quad\forall\phi$		(9)
	$\displaystyle\int_{\Omega}\frac{\partial\mathbf{u}}{\partial t}\cdot\bm{\phi}\;\mathrm{d}x+\mathcal{S}_{adv}+\mathcal{S}_{cor}+\mathcal{S}_{pg}+\mathcal{S}_{pa}$	$\displaystyle=$	$\displaystyle\mathcal{S}_{visc}+\mathcal{S}_{\bm{\tau}},\quad\forall\bm{\phi}.$		(10)

The bilinear forms of each term, $\mathcal{S}_{\cdot}$ , are defined below.

Let $\mathcal{T}$ stand for the triangulation of the domain $\Omega$ . The set of element interfaces is denoted by $\mathcal{I}=\{k\cap k^{\prime}|k,k^{\prime}\in\mathcal{T}\}$ , and $\mathbf{n}=(n_{x},n_{y})$ denotes the outward unit normal vector of an interface $e\in\mathcal{I}$ . On the interfaces, the DG₁ functions are discontinuous and do not have a unique value. We define the average and jump operators,

$\displaystyle\{\!\!\{a\}\!\!\}$	$\displaystyle=$	$\displaystyle\frac{1}{2}(a^{+}+a^{-}),$	(11)
$\displaystyle{[}\!{[}a\textbf{n}{]}\!{]}$	$\displaystyle=$	$\displaystyle a^{+}\mathbf{n}^{+}+a^{-}\mathbf{n}^{-},$	(12)
$\displaystyle{[}\!{[}\mathbf{u}\cdot\mathbf{n}{]}\!{]}$	$\displaystyle=$	$\displaystyle\mathbf{u}^{+}\cdot\mathbf{n}^{+}+\mathbf{u}^{-}\cdot\mathbf{n}^{-},$	(13)
$\displaystyle{[}\!{[}\mathbf{u}\mathbf{n}{]}\!{]}$	$\displaystyle=$	$\displaystyle\mathbf{u}^{+}\mathbf{n}^{+}+\mathbf{u}^{-}\mathbf{n}^{-},$	(14)

where the superscripts ‘ $+$ ’ and ‘ $-$ ’ arbitrarily label the values on either side of the interface, and $\mathbf{n}^{-}=-\mathbf{n}^{+}$ .

The $H\mathbf{u}$ divergence term is integrated by parts:

\displaystyle\mathcal{S}_{div}=-\int_{\Omega}H(\textbf{u}\cdot\nabla\phi)\;\mathrm{d}x+\int_{\mathcal{I}}(H^{*}\textbf{u}^{*})\cdot{[}\!{[}\phi\textbf{n}{]}\!{]}\;\mathrm{d}S+\int_{\Gamma_{o}}(H\textbf{u})\cdot\phi\textbf{n}\;\mathrm{d}S

(15)

where $H^{*}$ and $\textbf{u}^{*}$ denote the interface terms obtained from an approximate Riemann solver defined below. Note that the $H\mathbf{u}$ flux vanishes on the closed boundary $\Gamma_{c}$ .

The advection term is also integrated by parts (also here the flux across $\Gamma_{c}$ vanishes):

\displaystyle\mathcal{S}_{adv}=-\int_{\Omega}\nabla_{h}\cdot(\textbf{u}\bm{\phi})\cdot\textbf{u}\;\mathrm{d}x+\int_{\mathcal{I}}\{\!\!\{\textbf{u}\}\!\!\}\cdot{[}\!{[}\bm{\phi}(\textbf{u}\cdot\textbf{n}){]}\!{]}\;\mathrm{d}S+\int_{\Gamma_{o}}\textbf{u}\cdot\bm{\phi}(\textbf{u}\cdot\textbf{n})\;\mathrm{d}S.

(16)

The Coriolis and atmospheric pressure gradient terms read

	$\displaystyle\mathcal{S}_{cor}$	$\displaystyle=$	$\displaystyle\int_{\Omega}f\textbf{e}_{z}\wedge\textbf{u}\cdot\bm{\phi}\;\mathrm{d}x,$		(17)
	$\displaystyle\mathcal{S}_{pa}$	$\displaystyle=$	$\displaystyle\int_{\Omega}\frac{1}{\rho_{0}}\nabla p_{a}\cdot\bm{\phi}\;\mathrm{d}x.$		(18)

The pressure gradient term is integrated by parts

\displaystyle\mathcal{S}_{pg}=-\int_{\Omega}g\eta\nabla\cdot\bm{\phi}\;\mathrm{d}x+\int_{\mathcal{I}}g\eta^{*}{[}\!{[}\bm{\phi}\cdot\textbf{n}{]}\!{]}\;\mathrm{d}S+\int_{\Gamma}g\eta\bm{\phi}\cdot\textbf{n}\;\mathrm{d}S.

(19)

The viscosity operator is discretized with the interior penalty method [Rivière (\APACyear2008), Hillewaert (\APACyear2013)]:

$\displaystyle\mathcal{S}_{visc}$	$\displaystyle=$	$\displaystyle\int_{\Omega}(\nabla\bm{\psi}):\bm{\tau}_{\nu}\;\mathrm{d}x-\int_{\mathcal{I}}{[}\!{[}\bm{\psi}\textbf{n}{]}\!{]}\cdot\{\!\!\{\bm{\tau}_{\nu}\}\!\!\}\;\mathrm{d}S$	(20)
		$\displaystyle-\int_{\mathcal{I}}\{\!\!\{\nu_{h}\}\!\!\}{[}\!{[}\mathbf{u}\textbf{n}{]}\!{]}\cdot\{\!\!\{\nabla\bm{\psi}\}\!\!\}\;\mathrm{d}S+\int_{\mathcal{I}}\sigma\{\!\!\{\nu_{h}\}\!\!\}{[}\!{[}\mathbf{u}\textbf{n}{]}\!{]}\cdot{[}\!{[}\bm{\psi}\textbf{n}{]}\!{]}\;\mathrm{d}S$
		$\displaystyle-\int_{\Omega}\frac{\nu_{h}\nabla(H)}{H}\cdot(\nabla\mathbf{u}+(\nabla\mathbf{u})^{T})\cdot\bm{\phi}\;\mathrm{d}x$

where $\bm{\tau}_{\nu}=\nu_{h}\nabla\mathbf{u}$ denotes the viscous flux. The colon operator, $:$ , stands for the double inner product that contracts over both indices of the tensor expressions. The last term is an additional source term that accounts for the bathymetry gradient. The viscous flux is zero on $\Gamma_{c}$ and on $\Gamma_{o}$ .

Finally, the surface and bottom stress terms are given by

\displaystyle\mathcal{S}_{\bm{\tau}}=\int_{\Omega}\frac{\bm{\tau}_{w}+\bm{\tau}_{b}}{H\rho_{0}}\cdot\bm{\phi}\;\mathrm{d}x.

(21)

In this work, we use the Roe fluxes [Roe (\APACyear1981), LeVeque (\APACyear2002)] to stabilize the $\mathcal{S}_{div}$ and $\mathcal{S}_{pg}$ terms:

$\displaystyle\eta^{*}$	$\displaystyle=$	$\displaystyle\{\!\!\{\eta\}\!\!\}+\sqrt{\frac{g}{H}}{[}\!{[}\mathbf{u}\cdot\mathbf{n}{]}\!{]}$	(22)
$\displaystyle\mathbf{u}^{*}$	$\displaystyle=$	$\displaystyle\{\!\!\{\mathbf{u}\}\!\!\}+\sqrt{\frac{H}{g}}{[}\!{[}\eta\mathbf{n}{]}\!{]}$	(23)
$\displaystyle H^{*}$	$\displaystyle=$	$\displaystyle\eta^{*}+h.$	(24)

The spatial integrals are evaluated with a standard Gauss quadrature rule of degree 3. Denoting the individual DG₁ basis functions by $\phi_{i}$ and $\bm{\phi}_{i}$ , the discrete FE representation of $\eta$ and u fields are $\eta^{h}:=\sum_{i}\eta_{i}\phi_{i}$ and $\textbf{u}^{h}:=\sum_{i}\mathbf{u}_{i}\bm{\phi}_{i}$ , where $\eta_{i}$ and $\mathbf{u}_{i}$ are the nodal values. The superscript $h$ is used to denote the discrete FE fields in contrast to the continuous ones in (1-2). The model state (dual) vectors are $\bm{\eta}^{h}=(\eta_{1},\eta_{2},\ldots)$ , $\mathbf{u}^{h}=(u_{1},u_{2},\ldots,v_{1},v_{2},\ldots)$ ; the concatenated state vector is denoted by $\mathbf{q}^{h}=(\mathbf{u}^{h},\mathbf{\eta}^{h})$ . Replacing the test functions in (9-10) by $\phi_{i}$ and $\bm{\phi}_{i}$ , respectively, the equations can be written in a vector form:

\displaystyle\underline{\mathbf{M}}\frac{\partial\mathbf{q}^{h}}{\partial t}

\displaystyle=

\displaystyle\mathbf{S}(\mathbf{q}^{h}),

(25)

where $\underline{\mathbf{M}}$ denotes the 3-block DG₁ mass matrix, i.e. $\underline{\mathbf{M}}=\mathrm{diag}(\underline{\mathbf{M}}_{\phi},\underline{\mathbf{M}}_{\phi},\underline{\mathbf{M}}_{\phi})$ and $[\underline{\mathbf{M}}_{\phi}]_{i.j}=\int_{\Omega}\phi_{i}\phi_{j}dx$ , and $\mathbf{S}$ denotes the remaining bilinear forms. Hereafter we omit the $h$ superscript for brevity.

2.2 Time integration

The solution is marched in time with a fully implicit Runge-Kutta scheme. In this work, we use the two-stage, 2nd order accurate Diagonally Implicit Runge Kutta method DIRK(2,2), defined by the Butcher tableau [Ascher \BOthers. (\APACyear1997)],

\displaystyle\begin{array}[]{c|cc}c_{1}&a_{1,1}&0\\ c_{2}&a_{2,1}&a_{2,2}\\ \hline\cr&b_{1}&b_{2}\end{array}

\displaystyle\begin{array}[]{c|cc}\gamma&\gamma&0\\ 1&1-\gamma&\gamma\\ \hline\cr&1-\gamma&\gamma\end{array}

(32)

with $\gamma=(2-\sqrt{2})/2$ .

Denoting the time step by $\Delta t$ , solution at time $t^{n}$ by $\textbf{q}^{n}$ , and intermediate solutions by $\mathbf{q}^{(i)}$ , the $m$ -stage DIRK iteration reads

$\displaystyle\mathbf{q}^{(0)}$	$\displaystyle=$	$\displaystyle\mathbf{q}^{n},$	(33)
$\displaystyle\underline{\mathbf{M}}\mathbf{q}^{(i)}$	$\displaystyle=$	$\displaystyle\underline{\mathbf{M}}\mathbf{q}^{(i-1)}+\Delta t\sum_{j=1}^{i}a_{i,j}\mathbf{S}(\mathbf{q}^{(j)}),\quad\forall i=1,\ldots,m,$	(34)
$\displaystyle\mathbf{q}^{n+1}$	$\displaystyle=$	$\displaystyle\mathbf{q}^{(m)},$	(35)

For brevity, we denote the entire nonlinear forward model update operator by $F$ , i.e.

\displaystyle F(\mathbf{q}^{n+1},\mathbf{q}^{n})

\displaystyle=

\displaystyle 0.

(36)

Agglomerating the model state over all $N$ time steps into a single vector $\mathbf{Q}=(\mathbf{q}^{0},\ldots,\mathbf{q}^{N})$ , we can re-write the entire forward operator as

\displaystyle\mathcal{F}(\mathbf{Q})=0.

(37)

Note that the $\mathcal{F}$ operator is only introduced to simplify notation, it is never assembled. Indeed, $\mathcal{F}$ consists of nonlinear solves and it can only be evaluated by iterating over the time steps. As seen in (36), the forward update depends on past values of the model state. Let $\tilde{\mathcal{F}}$ denote a linearized forward operator. Due to the time dependency, $\tilde{\mathcal{F}}$ could be assembled into a lower-diagonal matrix operator. In what follows, we’ll see that the adjoint model requires a transpose of $\tilde{\mathcal{F}}$ , therefore reversing the time dependency.

3 Variational inverse problem

3.1 General description

The forward model depends on some unknown control parameters, $\bm{\theta}$ , that need to be estimated, i.e. $\mathcal{F}=\mathcal{F}(\mathbf{Q},\bm{\theta})$ . Let $J$ denote the user-defined, scalar-valued cost function we aim to minimize. The inverse optimization problem then reads

\displaystyle\min_{\mathbf{Q},\bm{\theta}}\,J(\mathbf{Q},\bm{\theta}),\quad\textrm{subject to}\,\,\mathcal{F}(\mathbf{Q},\bm{\theta})=0.

(38)

The minimization problem can be solved with standard gradient-based optimization methods provided that one can compute the gradient of $J$ . Following the notation of [Funke (\APACyear2013)] (where $\mathbf{Q}$ and $\bm{\theta}$ are defined as column vectors but the partial derivatives of a scalar, e.g., ${\partial J}/{\partial{\mathbf{Q}}}$ , are row vectors), the gradient can be written as:

\displaystyle\frac{\mathrm{d}J}{\mathrm{d}\bm{\theta}}=\frac{\partial J}{\partial\mathbf{Q}}\frac{\mathrm{d}\mathbf{Q}}{\mathrm{d}\bm{\theta}}+\frac{\partial J}{\partial\bm{\theta}}.

(39)

The latter term on the right hand side can be obtained by differentiating $J$ symbolically. The first term is unknown and can be computed from the discrete adjoint equation,

\displaystyle\left[\frac{\partial\tilde{\mathcal{F}}}{\partial\mathbf{Q}}\right]^{T}\mathbf{L}=\left[\frac{\partial J}{\partial\mathbf{Q}}\right]^{T},

(40)

where $\mathbf{L}=(\bm{\lambda}^{0},\bm{\lambda}^{1},\ldots,\bm{\lambda}^{N})$ is the time-agglomerated adjoint variable. The adjoint variables, $\bm{\lambda}^{n}$ , have the same dimension as the model state, $\textbf{q}^{n}$ . That is, for the model state variables, $\eta^{n}$ and $\textbf{u}^{n}$ , we have corresponding spatially and temporally varying adjoint fields. Note that due to the transpose, the adjoint equation must be solved backward in time (Fig. 2).

Once $\mathbf{L}$ is known, the gradient of $J$ can be computed:

\displaystyle\frac{\mathrm{d}J}{\mathrm{d}\bm{\theta}}

\displaystyle=

\displaystyle-\mathbf{L}^{T}\frac{\partial\tilde{\mathcal{F}}}{\partial\bm{\theta}}+\frac{\partial J}{\partial\bm{\theta}}.

(41)

The adjoint method is attractive as the cost of solving the adjoint equation, and hence its evaluation of the gradient, does not depend on the number of unknown parameters, $\bm{\theta}$ [Funke (\APACyear2013)]. This is in contrast to many other methods of calculating gradient information such as finite difference approximations, the tangent linear model, or ensemble methods, where computational cost increases with the number of control variables if evaluation of the full gradient is required.

With the pyadjoint library, the adjoint equation can be derived automatically [Farrell \BOthers. (\APACyear2013)]. Leveraging symbolic representation of the finite element discretization, (37), the adjoint equation is formed by differentiating the forward equations (Fig.2); a similar procedure is used to compute ${\partial\tilde{\mathcal{F}}}/{\partial\bm{\theta}}$ and ${\partial J}/{\partial\bm{\theta}}$ .

The time-dependent forward problem is solved first and pyadjoint is used to record the sequence of forward solve operations on “tape” (also known as the Wengert list, see e.g., \citeAgriewank2008) over the range of $N$ time steps. The adjoint equation (40) is solved backward in time by rewinding the tape and applying the linearized adjoint operators (Fig.2). In this work, the forward model state vectors, $\bm{Q}$ , are kept in memory and retrieved during the backward pass. For larger problems, this may exceed the amount of physical memory and $\bm{Q}$ must be stored to disk. Efficient pyadjoint checkpointing mechanism is under development. Once the adjoint variables, $\bm{\lambda}^{n}$ , are known, ${\mathrm{d}J}/{\mathrm{d}\bm{\theta}}$ can be evaluated.

Automated generation of the discrete adjoint model circumvents two major bottlenecks in adjoint modeling: First, implementing the adjoint model by hand is tedious, often comparable to the cost of implementing the forward model itself. Second, we utilize the same code generator as with the forward model to obtain low-level implementation for the adjoint model (Fig. 2). Thus the adjoint is computationally efficient, as the code generator can apply similar code optimization methods, such as loop unrolling or tiling, to the adjoint kernels.

Firedrake and the pyadjoint library provide an interface to derive and solve the adjoint model automatically. The user only needs to define the forward model $\mathcal{F}$ and the cost function $J$ . The gradient ${\partial J}/{\partial\bm{\theta}}$ can then be evaluated efficiently with only few modifications to the model code, enabling a wide range of inverse modeling applications.

Utilizing the gradient ${\mathrm{d}J}/{\mathrm{d}\bm{\theta}}$ , the optimization problem (38) is solved with a quasi-Newton method where the Hessian of $J$ is estimated during the iteration. The convergence of the quasi-Newton method depends on several factors, e.g. nonlinearity of the problem and smoothness of the cost function. In this work, we employ additional regularization to ensure smoothness of the spatially-variable control parameter field, $\bm{\theta}$ .

3.2 Water elevation optimization

In this work, we minimize the model misfit with respect to water elevation observations by varying the Manning bottom friction coefficient field.

Let $\eta_{o,i}^{n}$ , $i=1,\ldots,B$ , $n=1,\ldots,N$ denote the observation time series at $B$ tide gauges. The model elevation field is interpolated in space at the tide gauge locations, resulting in corresponding modeled time series $\eta_{m,i}^{n}$ . If a tide gauge lies outside the model domain, a nearest element center is used. Time-averaged time series are denoted by $\bar{\eta}_{o,i}=\sum_{n}\eta_{o,i}^{n}/N$ and the bias-removed (centered) time series by $\hat{\eta}_{o,i}^{n}=\eta_{o,i}^{n}-\bar{\eta}_{o,i}$ . The standard deviation of observations is given by $(\sigma_{o,i})^{2}=\sum_{n}(\hat{\eta}_{o,i}^{n})^{2}/N$ . The Centered Root Mean Square Deviation (CRMSD) is

\left(\mathrm{CRMSD}_{i}\right)^{2}=\frac{1}{N}\sum_{n=1}^{N}\left(\hat{\eta}_{m,i}^{n}-\hat{\eta}_{o,i}^{n}\right)^{2}.

(42)

The cost function is then defined as

J_{o}(\mathbf{Q},\bm{\theta})=\frac{1}{B}\sum_{i=1}^{B}\frac{1}{(\sigma_{o,i})^{2}}\left(\mathrm{CRMSD}_{i}\right)^{2}.

(43)

Note that the mistfit $J_{o}$ is computed with the centered time series $\hat{\eta}_{j,i}^{n},j=\{o,m\}$ . This is necessary because the model exhibits an SSH offset with respect to the observations that is not known a-priori. The model bias $\bar{\eta}_{m,i}-\bar{\eta}_{o,i}$ can be affected by the bathymetry (defined approximately with respect to the mean sea level), river discharge, the SSH bias at the open boundary, simulated water transport between basins, and the reference level of the observations. Similar SSH bias removal has used in previous data assimilation applications as well [Kurapov \BOthers. (\APACyear2011)]. Furthermore, the cost function is scaled by the inverse variance of the observations to ensure similar weight across the tidal and non-tidal tide gauges, for example. Indeed, the standard deviation of SSH exceeds 3 m in the English channel while only being a few centimeters in the Baltic Sea.

In this application, the control variable is the spatially varying Manning bottom friction coefficient field, i.e. $\bm{\theta}=\mu$ (see eq. (7)). The Manning field is discretized in space with continuous linear elements (Fig. 1 a). The optimization problem is ill-posed, i.e. there exist infinitely many Manning coefficient fields that minimize the cost function, $J_{o}$ (e.g., \citeAzhang2011). In addition, $J_{o}$ typically exhibits multiple local minima that prevent convergence to the global optimum. To this end, we use spatial regularization to penalize local variability of the $\bm{\theta}$ field, i.e. we augment the cost function with a regularization term $J_{r}(\bm{\theta})$ ,

	$\displaystyle J(\mathbf{Q},\bm{\theta})$	$\displaystyle=$	$\displaystyle J_{o}(\mathbf{Q},\bm{\theta})+J_{r}(\bm{\theta})$		(44)
	$\displaystyle J_{r}(\bm{\theta})$	$\displaystyle=$	$\displaystyle\alpha\\|\underline{\mathbf{H}}(\bm{\theta})(\Delta x)^{2}\\|_{2}^{2}=\alpha\int_{\Omega}\underline{\mathbf{H}}(\bm{\theta}):\underline{\mathbf{H}}(\bm{\theta})(\Delta x)^{4}\;\mathrm{d}x,$		(45)

where $\underline{\mathbf{H}}(\bm{\theta})$ denotes the Hessian of the control field with respect to the spatial coordinates, $\|\cdot\|_{2}$ is the $L^{2}$ norm, $\Delta x$ is the local mesh element size, and $\alpha$ is an unknown regularization parameter.

It should be noted that we only penalize the Hessian, i.e. the constant part and first gradients of $\bm{\theta}$ are not constrained at all. When computing $\underline{\mathbf{H}}(\bm{\theta})$ , we use a zero-gradient assumption for $\bm{\theta}$ at the boundary; the finite element implementation for computing $\underline{\mathbf{H}}(\bm{\theta})$ is presented in B. The penalty term is proportional to the Frobenius norm of the $\underline{\mathbf{H}}(\bm{\theta})$ (sum the squares of each element). As such, the regularization term is well-behaving in the sense that it is convex, positive definite, and invariant under coordinate system rotation.

The scaling by $\Delta x$ means that we effectively penalize the $\bm{\theta}$ variability within each element: $\partial\bm{\theta}/\partial x\Delta x\approx\Delta\bm{\theta}$ . This regularization term ensures that the Manning coefficient field remains smooth while higher variability is allowed in regions with high mesh resolution.

The optimization problem (38) is solved with the L-BFGS-B quasi-Newton method with bound constraints [Byrd \BOthers. (\APACyear1995)] from the SciPy package [Virtanen \BOthers. (\APACyear2020)]. During the iteration, we impose lower and upper bounds, $10^{-3}$ and $5^{-1}\ \mathrm{s}\ \mathrm{m}^{-1/3}$ , respectively, for the Manning coefficient field; in practice the upper bound is never reached.

To facilitate applications, an inverse modeling module was as added to Thetis. The module allows defining the cost function $J_{o}$ based on observational time series data at station locations and the regression term $J_{r}$ . The module also interfaces with Scipy’s L-BFGS-B method and provides file IO for the control and adjoint variables.

4 North Sea – Baltic Sea simulation

We apply the 2D shallow water model to the North Sea and Baltic Sea to simulate tidal and atmospherically-driven water elevation dynamics.

4.1 Model configuration

The model domain covers the North Sea and Baltic Sea (Fig. 3). The open boundary is placed beyond the continental shelf to allow reliable imposition of the tides due to greater water depth and weaker currents [Heemink \BOthers. (\APACyear2002), de Brye \BOthers. (\APACyear2010)]. In addition, the domain is sufficiently large so that most atmospherically-driven events, such as storm surges, can be generated within the model.

Fig. 4 presents the triangular unstructured mesh, generated with GMSH [Geuzaine \BBA Remacle (\APACyear2009)]. We first defined a scalar field indicating the desired mesh resolution across the domain. The coastal boundary line was extracted from the bathymetry raster at the 0.5 m contour. In order to generate a smooth coastline, the bathymetry was first smoothed with a Gaussian filter; stronger smoothing was applied in regions with coarser desired mesh resolution. This procedure ensures that the coastline is compatible with the bathymetric data and also appropriate for the intended mesh resolution. Higher mesh resolution is used along the coasts and especially in the Danish waters to better capture the complex geometry; resolution progressively decreases towards the open boundary.

Bathymetric data originates from the 1/16 arc minute EMODnet 2020 dataset [EMODnet Bathymetry Consortium (\APACyear2020)]. A second order P₂ discretization (Fig. 1 c) is used for the bathymetry field, $h$ , in order to improve the representation of small-scale features, such as narrow channels in the Danish Straits. During our initial tests we noted that the geometry and bathymetry of the channels play an important role in both the volume flux between the North Sea and Baltic Sea, as well as the reflection and refraction of tides in the region. The procedure of generating the P₂ bathymetry field is described in A. Wetting and drying is not used in the simulation as it has only a minor impact on the Baltic Sea SSH dynamics. Minimum water depth was set to 2 m. In areas with strong tides, the minimum depth was further increased to accomodate the tidal range (based on TPXO tidal model data).

At the open boundary, the model is forced with the TPXO 9 (v1) global tidal model [Egbert \BBA Erofeeva (\APACyear2002)] using all 15 constituents. We impose only tidal water elevation; water velocity is relaxed towards zero using the Roe fluxes. Our experiments indicated that including tidal velocity forcing does not improve model skill and may result in instabilities at the vicinity of the boundary. Subtidal SSH variation is not imposed. This configuration was found to be sufficient to represent both the mean SSH and atmospherically-driven SSH variability in the region of interest (the Danish Straits and the Baltic Sea).

A constant-in-time river discharge is imposed at 428 major rivers across the domain. The mean river discharge was computed from the EHYPE watershed model [Donnelly \BOthers. (\APACyear2015)] hindcast data.

Atmospheric wind stress and mean sea level pressure are imposed from the 2.5 km MetCoOp Ensemble Prediction System (MEPS) data obtained from the Norwegian Meteorological institute. The MEPS data set does not cover the western and southwestern parts of the model domain where the European Centre for Medium-range Weather Forecasts (ECMWF) HRES forecast data is used instead. The datasets are blended together in a 50 km overlapping region using a linear blending mask. Wind stress is computed with the [Large \BBA Yeager (\APACyear2008)] formulation from 10 m winds. The atmospheric pressure and wind stress fields are discretized with linear continuous P₁ elements (Fig. 1 a). Viscosity was set to a constant 5 $\textrm{m}^{2}/\textrm{s}$ throughout the domain.

The model time step was set to 1 h. The time step was chosen to minimize computational cost: A fully implicit solver allows longer time steps and can therefore reduce computational cost significantly [Kärnä \BOthers. (\APACyear2011)]. However, using a $\Delta t>1\ \textrm{h}$ resulted in slower convergence of the solver and thus higher overall cost. Using 1 h time step for modeling tidal dynamics is appropriate as we are using a second-order implicit Runge-Kutta solver which can represent nonlinear processes accurately. In addition, during the Runge-Kutta sub-iteration, all forcing fields are evaluated twice in each time step which, again, increases the accuracy. Our preliminary tests did not show any significant deterioration of the SSH performance with 1 h time step (not shown). We did notice, however, that the commonly-used, asymptotically unstable Crank-Nicolson time integrator did exhibit numerical instability manifested in a noisier velocity field.

The modeled SSH is compared against tide gauge observations at 56 tide gauges across the model domain (Fig. 3). The observational data was obtained from the Copernicus Marine Service (CMEMS) catalog [E.U. Copernicus Marine Service (\APACyear2015\APACexlab\BCnt2), E.U. Copernicus Marine Service (\APACyear2015\APACexlab\BCnt1)]. Only tide gauges that had nearly complete data coverage over the modeled period were included. The tide gauge time series were manually quality-checked to remove spurious SSH values (e.g., spikes).

4.2 Simulation period and initial condition

The Manning coefficient is optimized for a 16.5 day period, spanning from June 1 00:00 UTC to June 17 12:00 UTC, 2019. To exclude initial adjustment effects from the optimization, the cost function is not evaluated during the first 2.5 days, i.e. the model misfit is calculated over a 14-day period. This period is henceforth called the optimization period.

We chose to use a 14-day period to cover a full spring-neap tidal cycle as well as several seiche waves in the Baltic Sea. Early experiments showed that a shorter, 2-day period results in over-fitting as the model adjusts to particular realizations of coastal waves. Generally, choosing a suitable optimization window length is a trade-off, as the computational cost is directly proportional to the number of time steps.

On June 1, 2019, the model is initialized from a spun-up initial state. Using a realistic initial condition is important for the optimization process. According to our sensitivity runs, the mean SSH in the Baltic Sea reaches an equilibrium in roughly 1 month (not shown) and this transient adjustment must be excluded from the optimization process as it would skew the results. Moreover, a relatively good guess for the Manning field is needed for the spin-up run as bottom friction affects the volume flux to/from the Baltic Sea and therefore the mean SSH. The initial condition was generated as follows: First, 10 iterations of the optimization process are carried out starting the model from rest (zero initial water elevation and velocity; the Manning coefficient was set to 0.03 $\mathrm{s}\ \mathrm{m}^{-1/3}$ ). As a result, a first guess of a suitable Manning coefficient field, $\mu_{0}$ , is obtained. Using the $\mu_{0}$ field, a 1-month spin-up run (May 1, 2019 to June 1, 2019) is then carried out to generate an initial condition for June 1. The spin-up run used the ERA5 atmospheric reanalysis data as forcing [Hersbach \BOthers. (\APACyear2018)].

To initialize the final optimization, the Manning coefficient is set to 0.03 $\mathrm{s}\ \mathrm{m}^{-1/3}$ everywhere in the model domain. The optimization was run for 40 iterations, after which the changes in the friction field were small.

Once the optimal Manning coefficient field has been found, the results are validated with a 3-month simulation ranging from June 1 to August 31, 2019. The start date is the same as with the optimization period and the same model initial condition is used for elevation and velocity. The first 2.5 days are again excluded from the analysis. The longer validation run allows to verify that the Manning coefficient has not been tuned to fit particular events during the optimization period (over-fitting), i.e., the optimized model is able to represent SSH dynamics accurately in general.

4.3 Performance metrics

The model skill is quantified with standard statistical metrics. Using the notation from Section 3.2 (we drop the station index $i$ for brevity), the bias and correlation coefficient ( $R$ ) are defined as:

	$\displaystyle\mathrm{BIAS}$	$\displaystyle=$	$\displaystyle\bar{\eta}_{m}-\bar{\eta}_{o},$		(46)
	$\displaystyle R$	$\displaystyle=$	$\displaystyle\frac{1}{\sigma_{o}\sigma_{m}}\frac{1}{N}\sum_{n=1}^{N}\hat{\eta}_{m}^{n}\hat{\eta}_{o}^{n}.$		(47)

CRMSD is related to $\sigma_{m}$ and $R$ through the equation,

\displaystyle\mathrm{CRMSD}^{2}

\displaystyle=

\displaystyle\sigma_{o}^{2}+\sigma_{m}^{2}-\sigma_{o}\sigma_{m}R,

(48)

which can be visualized in a Taylor diagram [Taylor (\APACyear2001)]. In this work, we are normalizing the Taylor diagram by scaling the variables with $\sigma_{o}$ :

$\displaystyle\mathrm{NCRMSD}^{2}$	$\displaystyle=$	$\displaystyle 1+\sigma_{m}^{\prime 2}-\sigma_{m}^{\prime}R,$	(49)
$\displaystyle\mathrm{NCRMSD}$	$\displaystyle=$	$\displaystyle\frac{1}{\sigma_{o}}\mathrm{CRMSD},$	(50)
$\displaystyle\sigma_{m}^{\prime}$	$\displaystyle=$	$\displaystyle\frac{\sigma_{m}}{\sigma_{o}},$	(51)

where $\mathrm{NCRMSD}$ and $\sigma_{m}^{\prime}$ are the normalized CRMSD and standard deviation of the model, respectively. Normalization leads to dimensionless metrics and permits the comparison of different data sets.

5 Results

5.1 Optimization

The evolution of the cost function during the optimization is shown in Fig. 5. During the first 5 iterations the cost function decreases rapidly and begins to stagnate between iterations 20 and 40. The regularization parameter, $\alpha$ , was set to 400 as it appeared to result in a smooth Manning coefficient field without constraining the optimization too much.

The final optimized Manning coefficient field is shown in Fig. 6. The Manning field shows significant spatial variability throughout the model domain, i.e. not only in the vicinity of the tide gauge locations. This suggests that the optimization process alters the propagation of long surface waves instead of finding artificial local minima at the observation locations. Fig. 6 shows some oscillatory patterns in the Atlantic Ocean north of Scotland. These patterns are due to the spatial regularization method which does not ensure a completely smooth Manning coefficient field. In practice, however, the effect of these oscillations is small.

In the North Sea and around the British Isles, the friction coefficient is altered to improve the propagation and reflection of the tides; In the Danish waters (Fig. 6 b), friction is lowered in the Great Belt and across the Darss Sill. This is the dominant route for water transport between the Baltic Sea and the North Sea [Mohrholz \BOthers. (\APACyear2015), Gräwe \BOthers. (\APACyear2015), Mohrholz (\APACyear2018)] and the friction coefficient has a great impact on the barotropic volume flux, driven by water elevation difference across the Danish Straits.

In the Baltic Sea, the friction coefficient is close to the initial value especially in the deeper parts of the basin. In these regions (and also close to the open boundary) the gradient of the cost function shows very small values over the entire optimization process (not shown), suggesting that the cost function is not very sensitive to bottom friction in these regions. However, the Manning coefficient is altered at certain shallow regions, such as the Quark, Archipelago Sea and the Irbe Strait. These locations tend to have stronger depth-averaged currents and control the volume transport to the Bothnian Bay, Bothnian Sea, and Gulf of Riga, respectively. The Archipelago Sea is largely under-resolved with the used mesh resolution which plausibly explains the need for higher friction. \citeAmassmann2010 also found that adjoint-based sensitivity indicates higher bottom friction in the vicinity of unresolved islands. Overall, the optimized Manning coefficient field is complex and would be difficult to find through manual manipulation or brute force methods.

Example time series of the optimization process are shown in Fig. 7. In the Baltic Sea (panels a and b), the optimization mainly affects the slowly varying mean water elevation. This is related to the volumetric transport into the Baltic Sea. The Degerby station lies close to a nodal point of the basin’s seiche oscillations and thus it is often regarded as a metric of the total water volume in the basin. Panel (b) thus suggests that the emptying of the basin is initially underestimated; after 5 iterations, the model already tracks the observations fairly closely.

Panels (c) and (d) show examples of tidal stations in the Great Belt and Skagerrak, respectively. In both cases, the tidal range is initially overestimated but the model recovers a reasonable tidal range as the iteration progresses. The optimization also affects the phase of the tides to some extent. This is plausibly achieved by altering the relative strength of the tidal waves that propagate though the various parallel straits, for example. Based on our experience, these kinds of effects are difficult to achieve through manual manipulation of the Manning coefficient field. It should be noted that in the Danish Straits the solution converges slower compared to the Baltic stations (a and b) and the tidal waves are not fully replicated even after 40 iterations.

Figure 8 shows the estimated amplitude and phase of the M2 tidal constituent before and after the optimization (panels (a) and (b), respectively). Overall, the M2 tidal pattern is similar in both cases, but there are some significant differences. Before optimization, the model overestimates the M2 amplitude by several tens of centimeters in the west coast of Denmark and in the Danish Straits. Thus the overestimation of tidal range shown in Fig. 7 is a generic feature of the entire region. Optimization also affects the phase and moves the amphidromic points slightly. After optimization, the M2 amplitude and phase are in good agreement with the TPXO model (panel (c)).

5.2 Validation

Model performance metrics for the validation period are presented in Fig. 9. The blue and red markers indicate the initial and optimized model, respectively. The optimization improves the model performance at all stations. In the Baltic Sea, the CRMSD improves from 5 to 8 cm range to 2 to 6 cm range; In many locations in the Gotland Basin and the Bothnian Sea, the optimized model reaches CRMSD of 3 cm or less. Slightly higher CRMSD is observed in the Bothnian Bay and Gulf of Finland. In the Danish Straits, the deviation is higher, the optimized model’s CRMSD ranging between 4 and 8 cm. In this case, the improvement is substantial: At certain locations (e.g. Aarhus, Fredericia, Frederikshavn) the CRMSD is reduced by nearly a factor of 3. This is due to the complex reflection and interaction of tidal waves in this region which is difficult to replicate without careful tuning. For example, the tides propagate through the straits with different phase speeds, resulting in a complex superpositioning of tidal waves in the Great Belt; our initial tests indicate that a sufficiently good representation of the geometry is essential for capturing these processes.

The second panel in Fig. 9 shows the standard deviation of the observed and modeled and time series; The black markers indicate the observations. Also here, the optimization clearly improves the performance. The model has a slight tendency to underestimate variability at certain stations in the Gulf of Bothnia and Danish Straits.

Fig. 10 presents a Taylor diagram for the validation data. Also here the improvement in model skill is clear. In the Baltic Sea (round markers) the optimized model performs very well; at most stations the correlation coefficient is above 0.95 and standard deviation is close to the observations. In the Danish waters (star markers), the correlation coefficient is lower but still above 0.83 in all cases.

Two time series examples are shown in Fig. 11. Panel (a) shows water elevation at Degerby station for the entire validation period. Similarly to Fig. 7, the uncalibrated model tends to underestimate the slowly-varying SSH changes, suggesting that the transport through the Danish Straits is underestimated. The optimized model, on the other hand, tracks the volumetric changes quite accurately. Panel (b) shows tidal SSH signal at Korsor station in the Great Belt for a 7-day period outside the optimization window. Also in this case, the optimized model produces both a realistic tidal range as well as sub-tidal variability.

5.3 Computational cost

The simulations were carried out on AMD Rome 7H12 CPUs on the CSC Mahti cluster using 8 MPI processes. The forward model runs roughly 1800 times faster than real time: a 1-day and a 1-month simulation take roughly 48 s and 24 min, respectively. Evaluating the gradient of $J$ (including both the forward and adjoint solves) in the 16.5 day optimization period takes roughly 45 min, i.e. about 3.4 times longer than running the forward model alone. In each L-BFGS-B iteration typically only one, and at most two, gradient evaluations are needed. Running the entire optimization for 40 iterations took 36 h.

Comparing the cost against other adjoint model implementations is not straightforward as computational cost is rarely reported. \citeAyaremchuk2016 used a manually-coded adjoint of the Navy Coastal Ocean Model (NCOM); a single tangent linear model (TLM) and adjoint evaluation costs approximately 11 times as much as the forward model. \citeAvillaret2016 used an AD-generated TLM to study the sensitivity of the Telemac2D/Sisyphe morphodynamic model. They report that one TLM evaluation costs about 3 times as much as the forward model. Considering that a TLM evaluation is cheaper than the full forward-backward solve, this suggests that our adjoint implementation is at least as efficient as the implementations referenced above.

6 Discussion

We have presented an adjoint-based optimization of a 2D water elevation model. One of our main findings is that for the optimization process to work, the model misfit must be dominated by the control variable, i.e. bottom friction in this case. Otherwise, the optimization may end up compensating errors due to other sources, e.g., by generating artificial friction patterns near the boundary in the case of poor boundary forcings. Consequently, one needs a sufficiently good mesh, bathymetry, boundary conditions and atmospheric forcing. The open boundaries should be located sufficiently far away from the domain of interest. In addition, a realistic initial condition is crucial for robust optimization.

In this work we are using a 2D shallow water model. As such, our model does not include water density effects, or baroclinic processes. Wind wave effects are also not included. Water density difference between the North Sea and the northern Baltic Sea introduces a steady elevation gradient across the Baltic basin. Baroclinic processes can be important especially in the Danish Straits where density gradients are strong and fluctuate with in- and outflow events. The freshwater outflow in Kattegat/Skagerrak can exhibit baroclinic eddies, which affect local water elevation. Wind wave effects are most pronounced under storm conditions and can enhance storm surges via Stokes drift or wave build-up near the coast [Kanarik \BOthers. (\APACyear2021)]. The results indicate, however, that a 2D barotropic model can replicate the majority of the SSH variability with good accuracy.

The presented inverse modeling procedure resembles the so-called four-dimensional variational assimilation (4D-Var) method widely used in both atmospheric [Le Dimet \BBA Talagrand (\APACyear1986)] and oceanic [Thacker \BBA Long (\APACyear1988), Wunsch \BBA Heimbach (\APACyear2007)] applications. Compared to simpler methods (e.g. 3D-Var), the time dependency is consistently treated with the forward-backward solution procedure. In addition, the constraint (37) ensures that the solution always satisfies the model equations which is not the case in many widely-used data-assimilation methods (e.g. simple nudging of the model state, and filtering methods including the Kalman filter and its approximations). One notable difference is that some data-assimilation implementations use the ensemble method to estimate the gradient of the cost function, which is generally more expensive than solving the adjoint equation as the cost is directly proportional to the size of the ensemble.

One common pitfall in gradient-based optimization is over-fitting which must be controlled with some form of regularization. \citeAheemink2002 and \citeAzhang2011 reduced the dimensionality of the friction field considering only a small subset of the grid’s nodal values. In this work, in addition to using a sufficiently long optimization period, we included a regularization term to penalize the second derivative of the control variable. This procedure is sufficient to avoid strong “bipolar” friction adjustment at the tide gauge locations. The final Manning coefficient field is complex and highly variable across the domain; friction is altered also in regions where no observation data is used (e.g., around the British Isles and the Irish Sea). It is worth noting, however, that shallow water surface waves are typically relatively long and smooth which also reduces the risk of spatial over-fitting.

The presented optimization process is robust as most tests resulted in a similar friction field. We did observe some dependency on the mesh, however. As such, the optimized friction field should be regarded as a model and mesh dependent parameter with limited physical interpretability (e.g. with respect to sea bed roughness). Indeed, capturing SSH dynamics requires suitable dissipation at the right locations and the model’s dissipation characteristics are governed by the discretization, mesh resolution as well as physical parametrizations.

The fact that the model skill is lower in the Danish waters compared to the Baltic Sea stations is consistent with other modeling results. Compared to the 1 nautical mile Nemo-Nordic 2.0 3D baroclinic model [Kärnä \BOthers. (\APACyear2021)], the presented results are comparable or better. Although the simulation period, model configuration, and forcings are different, \citeAkarna2021 state similar model skill in the Baltic Sea and Danish waters (CRMSD below 7 cm and 10 cm, respectively). However, at certain locations, e.g. Fredericia, accuracy is lower, CRMSD being above 12 cm. This is most likely due to the fact that the complex geometry of the Danish Straits cannot be properly represented with a 1 nautical mile (1.8 km) structured grid.

As tuning of bottom friction is necessary in many coastal applications, the presented optimization procedure is a promising step towards automated calibration of ocean models. Manual tuning of spatially variable friction coefficient can be a very time consuming task. Without an adjoint model, the options for automated tuning are limited: one has to rely on a brute-force search, or ensemble-based estimates of the cost function gradient, for example. In both approaches, computational cost is proportional to the degrees of freedom (number of nodes) of the control field. In this application, one iteration takes 3.4 times longer than running the forward model, i.e. the cost is roughly equivalent to running a 3-member ensemble. The accuracy, however, is much better as the gradient ${\mathrm{d}J}/{\mathrm{d}\bm{\theta}}$ is computed exactly.

Adjoint-based optimization has many similarities with machine learning methods [Sonnewald \BOthers. (\APACyear2021)]. The backward-in-time adjoint iteration resembles the backpropagation training methods. In addition, regularization is a key ingredient to avoid over-fitting and a separate validation set is commonly used to inform the final model selection. In this application, the cost function is defined as the centralized root mean square deviation, scaled by the observation variance. This is in fact equivalent to whitening the target signal (i.e. removing the bias and scaling to unit variance), commonly used as a pre-processing step in machine learning. In contrast to data-driven machine learning methods, the benefit of adjoint-based physical modeling is that the results are guaranteed to be physically sound, i.e. satisfy the physical principles encoded in the forward equations.

7 Conclusions

We have presented an adjoint-based bottom friction optimization in a 2D water elevation model for the Baltic Sea. The discrete adjoint model is automatically generated from the symbolic forward model equations. The Manning coefficient optimization procedure is robust and yields good, generalizable results. Achieving similar performance via manual tuning would be challenging. The model performs well, especially in the complex Danish Straits region that is difficult to model with structured grid models (e.g., \citeAkarna2021), highlighting the benefit of variable-resolution unstructured meshes.

Leveraging the symbolic representation of the forward model, the adjoint model can be generated and solved with very few user efforts: The user only needs to implement the forward model and the desired cost function. The resulting adjoint model is computationally efficient: evaluating the gradient of the cost function takes 3.4 times as long as running the forward model alone. This work demonstrates that domain specific language frameworks, based on high-level abstractions and code generation, can facilitate automated model analysis and provide easy access to generic inverse modeling capability.

Appendix A Quadratic bathymetry field

As the coastal topography exerts a leading-order control on water circulation, it is crucial to obtain an accurate numerical representation of the bathymetry. Typically the accuracy of the bathymetric field is restricted by the mesh resolution. In this work, we discretize the bathymetry field with a second order P₂ element (Fig. 1 c) to increase the intrinsic resolution.

A common bottleneck with order $p>1$ polynomials is the loss of monotonicity near sharp gradients: The polynomial representation can lead to large oscillations that exceed the bounds of the original data. Water depth may become negative in the vicinity of steep slopes, for example.

We generate a smooth P₂ bathymetry field as follows: First, we generate a refined triangular mesh by splitting each triangle uniformly into four (i.e. joining the P₂ element nodes). On the refined mesh, we define a preliminary P₁ bathymetry field and interpolate the raw raster bathymetry data on it. This field does not exhibit any spurious overshoots due to the linear basis functions. The preliminary bathymetry field is then smoothed with Laplacian diffusion. Given an preliminary bathymetry field $h^{(1)}$ and a test function $\phi$ in the P₁ space, the smoothing operator is

	$\displaystyle r$	$\displaystyle=$	$\displaystyle\frac{(\|\nabla h^{(n)}\|\Delta x)^{\gamma}}{2(h^{(n)})^{\beta}}$		(52)
	$\displaystyle\int_{\Omega}(h^{(n+1)}-h^{(n)})\phi dx$	$\displaystyle=$	$\displaystyle-\int_{\Omega}\alpha r\nabla h^{(n)}\cdot\nabla\phi dx,\quad\forall\phi$		(53)

with $\alpha=20$ , $\gamma=5/2$ , and $\beta=1/2$ . On the right hand side, the Laplacian diffusion operator has been integrated by parts, omitting the interface terms. The diffusion coefficient, $r$ , depends on the local bathymetry change in each element. It resembles the dimensionless bathymetry slope metric commonly used in sigma coordinate models (recovered by setting $\gamma=\beta=1$ ; \citeAhaney1991,mellor1998).

The smoother is applied 30 times. We then use $L^{2}$ projection to cast the smoothed preliminary bathymetry on a P₂ field on the original mesh. Although the procedure is not guaranteed to be monotonic, the smoothing step essentially ensures that the projection does not generate severe overshoots. A comparison of the final P₂ bathymetry and bathymetry interpolated on the P₁ elements on the original mesh are shown in Fig. 12. The P₂ discretization allows more realistic representation of the narrow channels in the Danish Straits. In the P₁ case, the channels are not entirely continuous which affects volume transport.

Appendix B Hessian recovery procedure

The cost function described in Section 3.2 involves a regularization term that depends on the Hessian of the control field with respect to the spatial coordinates. In practice, the control field is discretized using linear finite elements, which are not twice continuously differentiable. In order to approximate the second spatial derivatives of the (continuous) control field, we apply an $L^{2}$ projection recovery procedure. The gradient of the discrete field is projected into $L^{2}$ space. A second application then obtains the Hessian approximation. We opt to perform the Hessian recovery in P₁ spaces and apply the two recovery steps simultaneously using a mixed finite element method.

Let $\theta^{h}$ denote the finite element approximation of the control field and $\mathbf{g}^{h}(\theta)$ and $\underline{\mathbf{H}}^{h}(\theta)$ denote the corresponding approximations of its gradient and Hessian, the latter of which we seek to obtain. Let $\bm{\phi}$ and $\underline{\bm{\phi}}$ denote test functions in vector and tensor P₁ spaces, respectively. The mixed formulation may then be written as

	$\displaystyle\int_{\Omega}\mathbf{g}^{h}\cdot\bm{\phi}\;\mathrm{d}x=-\int_{\Omega}\theta^{h}(\nabla\cdot\bm{\phi})\;\mathrm{d}x+\int_{\Gamma}\theta^{h}(\bm{\phi}\cdot\mathbf{n})\;\mathrm{d}s+\int_{\mathcal{I}}\{\!\!\{\theta^{h}\}\!\!\}{[}\!{[}\bm{\phi}\cdot\mathbf{n}{]}\!{]}\;\mathrm{d}S\quad\forall\bm{\phi}$		(54)
	$\displaystyle\int_{\Omega}\underline{\mathbf{H}}^{h}:\underline{\bm{\phi}}\;\mathrm{d}x=-\int_{\Omega}\mathbf{g}^{h}\cdot(\nabla\cdot\underline{\bm{\phi}})\;\mathrm{d}x+\int_{\Gamma}\mathbf{g}^{h}\cdot(\underline{\bm{\phi}}\cdot\mathbf{n})\;\mathrm{d}s+\int_{\mathcal{I}}\{\!\!\{\mathbf{g}^{h}\}\!\!\}\cdot{[}\!{[}\underline{\bm{\phi}}\cdot\mathbf{n}{]}\!{]}\;\mathrm{d}S\quad\forall\underline{\bm{\phi}}$		(55)

where $\int_{\Omega}\underline{\mathbf{A}}:\underline{\mathbf{B}}\;\mathrm{d}x$ denotes the standard $L^{2}$ inner product on $\Omega$ for two tensor-valued functions $\underline{\mathbf{A}}$ and $\underline{\mathbf{B}}$ . An advantage of the formulation in (54)–(55) is that it does not require the control field to have any spatial derivatives at all, meaning it can be applied to both continuous and discontinuous fields.

Taken together, (54) and (55) give rise to a $2\times 2$ block matrix system. We solve this system using Schur complement preconditioners and GMRES as the linear solver. A major advantage of our approach is that the cost function $J$ is defined symbolically, including the computation of the Hessian $\underline{\mathbf{H}}(\theta)$ . As such, the solution of (54)–(55) to approximate the Hessian is automatically taken into account in the computation of the gradient $\mathrm{d}J/\mathrm{d}\bm{\theta}$ during the optimization process.

Open Research

The model source code is publicly available. Firedrake, and its components, may be obtained from https://www.firedrakeproject.org (last access: 22 June 2023); Thetis may be obtained from http://thetisproject.org (last access: 22 June 2023). The exact software versions used to produce the results in this paper have been archived on Zenodo [zenodo/Firedrake (\APACyear2022), zenodo/Thetis (\APACyear2022)]. The model configuration as well as the model and observation time series used to produce the presented results are also available on Zenodo [zenodo/Data (\APACyear2022)]. Tidal forcing and inversion were computed with the Uptide Python package [zenodo/Uptide (\APACyear2020)]. Data analysis and visualization were carried out with Matplotlib [Hunter (\APACyear2007)] and Iris [Met Office (\APACyear2010 - 2022)] Python packages.

The MEPS atmospheric forecast data can be obtained from the Norwegian Meteorological Institute, https://github.com/metno/NWPdocs/wiki/Data-access (last access: 22 June 2023). The ECMWF (European Centre for Medium-Range Weather Forecasts) HRES atmospheric forecast data can be accessed from the ECMWF services, subject to license restrictions, https://www.ecmwf.int/en/forecasts/accessing-forecasts (last access: 22 June 2023). The ERA5 reanalysis atmospheric data can be downloaded from the Copernicus Climate Change Service (C3S) Climate Data Store (CDS) [Hersbach \BOthers. (\APACyear2018)]. The EHYPE river runoff data can be obtained from the Swedish Meteorological and Hydrological Institute, subject to license restrictions. The TPXO global tidal atlas (TPXO9-atlas-v1; \citeAegbert2002) can be obtained from https://www.tpxo.net (last access: 22 June 2023). Bathymetric data is provided by the \citeAemodnet2020. The observational data can be obtained from the Copernicus Marine Service [E.U. Copernicus Marine Service (\APACyear2015\APACexlab\BCnt2), E.U. Copernicus Marine Service (\APACyear2015\APACexlab\BCnt1)].

Acknowledgements.

The authors wish to acknowledge CSC – IT Center for Science, Finland, for computational resources. The work has received support from the Project HPC-EUROPA3 (INFRAIA-2016-1-730897), with the support of the EC Research Innovation Action under the H2020 Programme; in particular, J. G. Wallwork gratefully acknowledges the support of the Finnish Meteorological Institute and the computer resources and technical support provided by CSC. J. G. Wallwork was funded under the embedded CSE programme of the ARCHER2 UK National Supercomputing Service (http://www.archer2.ac.uk). S. C. Kramer acknowledges support from EPSRC under grant EP/R029423/1. We are grateful to Prof. Matthew Piggott (Imperial College London) for his insightful comments and suggestions. We thank the four anonymous reviewers for their comments that helped improve the manuscript.

References

Almeida \BOthers. (\APACyear2018) \APACinsertmetastaralmeida2018{APACrefauthors}Almeida, T\BPBIG., Walker, D\BPBIT.\BCBL \BBA Warnock, A\BPBIM. \APACrefYearMonthDay2018. \BBOQ\APACrefatitleEstimating River Bathymetry from Surface Velocity Observations Using Variational Inverse Modeling Estimating river bathymetry from surface velocity observations using variational inverse modeling.\BBCQ \APACjournalVolNumPagesJournal of Atmospheric and Oceanic Technology35121–34. {APACrefDOI} 10.1175/jtech-d-17-0075.1 \PrintBackRefs\CurrentBib
Alnæs \BOthers. (\APACyear2014) \APACinsertmetastaralnaes2014{APACrefauthors}Alnæs, M\BPBIS., Logg, A., Ølgaard, K\BPBIB., Rognes, M\BPBIE.\BCBL \BBA Wells, G\BPBIN. \APACrefYearMonthDay2014. \BBOQ\APACrefatitleUnified form language Unified form language.\BBCQ \APACjournalVolNumPagesACM Transactions on Mathematical Software4021–37. {APACrefDOI} 10.1145/2566630 \PrintBackRefs\CurrentBib
Ascher \BOthers. (\APACyear1997) \APACinsertmetastarascher1997{APACrefauthors}Ascher, U\BPBIM., Ruuth, S\BPBIJ.\BCBL \BBA Spiteri, R\BPBIJ. \APACrefYearMonthDay1997. \BBOQ\APACrefatitleImplicit-explicit Runge-Kutta methods for time-dependent partial differential equations Implicit-explicit runge-kutta methods for time-dependent partial differential equations.\BBCQ \APACjournalVolNumPagesApplied Numerical Mathematics252-3151–167. {APACrefDOI} 10.1016/s0168-9274(97)00056-1 \PrintBackRefs\CurrentBib
Balay \BOthers. (\APACyear2021) \APACinsertmetastarpetsc2021{APACrefauthors}Balay, S., Abhyankar, S., Adams, M\BPBIF., Benson, S., Brown, J., Brune, P.\BDBLZhang, J. \APACrefYearMonthDay2021. \APACrefbtitlePETSc/TAO Users Manual PETSc/TAO users manual \APACbVolEdTR\BTR \BNUM ANL-21/39 - Revision 3.16. \APACaddressInstitutionArgonne National Laboratory. \PrintBackRefs\CurrentBib
Byrd \BOthers. (\APACyear1995) \APACinsertmetastarbyrd1995{APACrefauthors}Byrd, R\BPBIH., Lu, P., Nocedal, J.\BCBL \BBA Zhu, C. \APACrefYearMonthDay1995. \BBOQ\APACrefatitleA Limited Memory Algorithm for Bound Constrained Optimization A limited memory algorithm for bound constrained optimization.\BBCQ \APACjournalVolNumPagesSIAM Journal on Scientific Computing1651190–1208. {APACrefDOI} 10.1137/0916069 \PrintBackRefs\CurrentBib
Clare \BOthers. (\APACyear2022) \APACinsertmetastarclare2022{APACrefauthors}Clare, M\BPBIC., Kramer, S\BPBIC., Cotter, C\BPBIJ.\BCBL \BBA Piggott, M\BPBID. \APACrefYearMonthDay2022. \BBOQ\APACrefatitleCalibration, inversion and sensitivity analysis for hydro-morphodynamic models Calibration, inversion and sensitivity analysis for hydro-morphodynamic models.\BBCQ \APACjournalVolNumPagesComputers and Geosciencesjust accepted - should have a complete ref shortly. \PrintBackRefs\CurrentBib
Daewel \BBA Schrum (\APACyear2013) \APACinsertmetastardaewel2013{APACrefauthors}Daewel, U.\BCBT \BBA Schrum, C. \APACrefYearMonthDay2013. \BBOQ\APACrefatitleSimulating long-term dynamics of the coupled North Sea and Baltic Sea ecosystem with ECOSMO II: Model description and validation Simulating long-term dynamics of the coupled north sea and baltic sea ecosystem with ECOSMO II: Model description and validation.\BBCQ \APACjournalVolNumPagesJournal of Marine Systems119-12030–49. {APACrefDOI} 10.1016/j.jmarsys.2013.03.008 \PrintBackRefs\CurrentBib
de Brye \BOthers. (\APACyear2010) \APACinsertmetastardebrye2010{APACrefauthors}de Brye, B., de Brauwere, A., Gourgue, O., Kärnä, T., Lambrechts, J., Comblen, R.\BCBL \BBA Deleersnijder, E. \APACrefYearMonthDay2010. \BBOQ\APACrefatitleA finite-element, multi-scale model of the Scheldt tributaries, river, estuary and ROFI A finite-element, multi-scale model of the Scheldt tributaries, river, estuary and ROFI.\BBCQ \APACjournalVolNumPagesCoastal Engineering579850–863. {APACrefDOI} 10.1016/j.coastaleng.2010.04.001 \PrintBackRefs\CurrentBib
Donnelly \BOthers. (\APACyear2015) \APACinsertmetastardonnelly2015{APACrefauthors}Donnelly, C., Andersson, J\BPBIC.\BCBL \BBA Arheimer, B. \APACrefYearMonthDay2015. \BBOQ\APACrefatitleUsing flow signatures and catchment similarities to evaluate the E-HYPE multi-basin model across Europe Using flow signatures and catchment similarities to evaluate the e-HYPE multi-basin model across europe.\BBCQ \APACjournalVolNumPagesHydrological Sciences Journal612255–273. {APACrefDOI} 10.1080/02626667.2015.1027710 \PrintBackRefs\CurrentBib
Egbert \BBA Erofeeva (\APACyear2002) \APACinsertmetastaregbert2002{APACrefauthors}Egbert, G\BPBID.\BCBT \BBA Erofeeva, S\BPBIY. \APACrefYearMonthDay2002. \BBOQ\APACrefatitleEfficient Inverse Modeling of Barotropic Ocean Tides Efficient inverse modeling of barotropic ocean tides.\BBCQ \APACjournalVolNumPagesJournal of Atmospheric and Oceanic Technology192183–204. {APACrefDOI} 10.1175/1520-0426(2002)019¡0183:eimobo¿2.0.co;2 \PrintBackRefs\CurrentBib
EMODnet Bathymetry Consortium (\APACyear2020) \APACinsertmetastaremodnet2020{APACrefauthors}EMODnet Bathymetry Consortium. \APACrefYearMonthDay2020. \APACrefbtitleEMODnet Digital Bathymetry (DTM 2020). Emodnet digital bathymetry (dtm 2020). \APACaddressPublisherEuropean Marine Observation and Data Network (EMODnet). {APACrefDOI} https://doi.org/10.12770/bb6a87dd-e579-4036-abe1-e649cea9881a \PrintBackRefs\CurrentBib
E.U. Copernicus Marine Service (\APACyear2015\APACexlab\BCnt1) \APACinsertmetastarcmems-northsea{APACrefauthors}E.U. Copernicus Marine Service. \APACrefYearMonthDay2015\BCnt1. \APACrefbtitleAtlantic- European North West Shelf - Ocean In-Situ Near Real Time observations. Atlantic- European North West Shelf - ocean in-situ near real time observations. \APACaddressPublisherE.U. Copernicus Marine Service. {APACrefDOI} 10.48670/MOI-00045 \PrintBackRefs\CurrentBib
E.U. Copernicus Marine Service (\APACyear2015\APACexlab\BCnt2) \APACinsertmetastarcmems-baltic{APACrefauthors}E.U. Copernicus Marine Service. \APACrefYearMonthDay2015\BCnt2. \APACrefbtitleBaltic Sea - In Situ Near Real Time Observations. Baltic Sea - in situ near real time observations. \APACaddressPublisherE.U. Copernicus Marine Service. {APACrefDOI} 10.48670/MOI-00032 \PrintBackRefs\CurrentBib
Farrell \BOthers. (\APACyear2013) \APACinsertmetastarfarrell2013{APACrefauthors}Farrell, P\BPBIE., Ham, D\BPBIA., Funke, S\BPBIW.\BCBL \BBA Rognes, M\BPBIE. \APACrefYearMonthDay2013. \BBOQ\APACrefatitleAutomated Derivation of the Adjoint of High-Level Transient Finite Element Programs Automated derivation of the adjoint of high-level transient finite element programs.\BBCQ \APACjournalVolNumPagesSIAM Journal on Scientific Computing354C369–C393. {APACrefDOI} 10.1137/120873558 \PrintBackRefs\CurrentBib
Funke (\APACyear2013) \APACinsertmetastarfunke2013{APACrefauthors}Funke, S. \APACrefYear2013. \APACrefbtitleThe automation of PDE-constrained optimisation and its applications The automation of PDE-constrained optimisation and its applications \APACtypeAddressSchool\BPhDImperial College London London. {APACrefDOI} 10.25560/11079 \PrintBackRefs\CurrentBib
Funke \BOthers. (\APACyear2014) \APACinsertmetastarfunke2014{APACrefauthors}Funke, S\BPBIW., Farrell, P\BPBIE.\BCBL \BBA Piggott, M. \APACrefYearMonthDay2014. \BBOQ\APACrefatitleTidal turbine array optimisation using the adjoint approach Tidal turbine array optimisation using the adjoint approach.\BBCQ \APACjournalVolNumPagesRenewable Energy63658–673. \PrintBackRefs\CurrentBib
Funke \BOthers. (\APACyear2016) \APACinsertmetastarfunke2016{APACrefauthors}Funke, S\BPBIW., Kramer, S\BPBIC.\BCBL \BBA Piggott, M\BPBID. \APACrefYearMonthDay2016. \BBOQ\APACrefatitleDesign optimisation and resource assessment for tidal-stream renewable energy farms using a new continuous turbine approach Design optimisation and resource assessment for tidal-stream renewable energy farms using a new continuous turbine approach.\BBCQ \APACjournalVolNumPagesRenewable energy991046–1061. \PrintBackRefs\CurrentBib
Geuzaine \BBA Remacle (\APACyear2009) \APACinsertmetastargeuzaine2009{APACrefauthors}Geuzaine, C.\BCBT \BBA Remacle, J\BHBIF. \APACrefYearMonthDay2009. \BBOQ\APACrefatitleGmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities.\BBCQ \APACjournalVolNumPagesInternational Journal for Numerical Methods in Engineering79111309–1331. {APACrefDOI} 10.1002/nme.2579 \PrintBackRefs\CurrentBib
Goss \BOthers. (\APACyear2021) \APACinsertmetastargoss2021{APACrefauthors}Goss, Z., Coles, D., Kramer, S.\BCBL \BBA Piggott, M. \APACrefYearMonthDay2021. \BBOQ\APACrefatitleEfficient economic optimisation of large-scale tidal stream arrays Efficient economic optimisation of large-scale tidal stream arrays.\BBCQ \APACjournalVolNumPagesApplied Energy295116975. \PrintBackRefs\CurrentBib
Gräwe \BOthers. (\APACyear2015) \APACinsertmetastargrawe2015{APACrefauthors}Gräwe, U., Naumann, M., Mohrholz, V.\BCBL \BBA Burchard, H. \APACrefYearMonthDay2015. \BBOQ\APACrefatitleAnatomizing one of the largest saltwater inflows into the Baltic Sea in December 2014 Anatomizing one of the largest saltwater inflows into the Baltic Sea in December 2014.\BBCQ \APACjournalVolNumPagesJournal of Geophysical Research: Oceans120117676–7697. {APACrefDOI} 10.1002/2015jc011269 \PrintBackRefs\CurrentBib
Griewank \BBA Walther (\APACyear2008) \APACinsertmetastargriewank2008{APACrefauthors}Griewank, A.\BCBT \BBA Walther, A. \APACrefYear2008. \APACrefbtitleEvaluating derivatives: principles and techniques of algorithmic differentiation Evaluating derivatives: principles and techniques of algorithmic differentiation. \APACaddressPublisherSIAM. \PrintBackRefs\CurrentBib
Haney (\APACyear1991) \APACinsertmetastarhaney1991{APACrefauthors}Haney, R\BPBIL. \APACrefYearMonthDay1991. \BBOQ\APACrefatitleOn the Pressure Gradient Force over Steep Topography in Sigma Coordinate Ocean Models On the pressure gradient force over steep topography in sigma coordinate ocean models.\BBCQ \APACjournalVolNumPagesJournal of Physical Oceanography214610–619. {APACrefDOI} 10.1175/1520-0485(1991)021¡0610:otpgfo¿2.0.co;2 \PrintBackRefs\CurrentBib
Heemink \BOthers. (\APACyear2002) \APACinsertmetastarheemink2002{APACrefauthors}Heemink, A., Mouthaan, E., Roest, M., Vollebregt, E., Robaczewska, K.\BCBL \BBA Verlaan, M. \APACrefYearMonthDay2002. \BBOQ\APACrefatitleInverse 3D shallow water flow modelling of the continental shelf Inverse 3D shallow water flow modelling of the continental shelf.\BBCQ \APACjournalVolNumPagesContinental Shelf Research223465–484. {APACrefDOI} 10.1016/s0278-4343(01)00071-1 \PrintBackRefs\CurrentBib
Hersbach \BOthers. (\APACyear2018) \APACinsertmetastarera5{APACrefauthors}Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horányi, A., Muñoz Sabater, J.\BDBLThépaut, J\BHBIN. \APACrefYearMonthDay2018. \APACrefbtitleERA5 hourly data on single levels from 1979 to present. ERA5 hourly data on single levels from 1979 to present. \APACaddressPublisherCopernicus Climate Change Service (C3S) Climate Data Store (CDS). {APACrefDOI} 10.24381/cds.adbb2d47 \PrintBackRefs\CurrentBib
Hillewaert (\APACyear2013) \APACinsertmetastarhillewaert2013{APACrefauthors}Hillewaert, K. \APACrefYear2013. \APACrefbtitleDevelopment of the discontinuous Galerkin method for high-resolution, large scale CFD and acoustics in industrial geometries Development of the discontinuous Galerkin method for high-resolution, large scale CFD and acoustics in industrial geometries \APACtypeAddressSchool\BUPhD. \APACaddressSchoolUniversité catholique de Louvain. \PrintBackRefs\CurrentBib
Homolya \BOthers. (\APACyear2018) \APACinsertmetastarhomolya2018{APACrefauthors}Homolya, M., Mitchell, L., Luporini, F.\BCBL \BBA Ham, D\BPBIA. \APACrefYearMonthDay2018. \BBOQ\APACrefatitleTSFC: A Structure-Preserving Form Compiler TSFC: A structure-preserving form compiler.\BBCQ \APACjournalVolNumPagesSIAM Journal on Scientific Computing403C401–C428. {APACrefDOI} 10.1137/17m1130642 \PrintBackRefs\CurrentBib
Hordoir \BOthers. (\APACyear2019) \APACinsertmetastarhordoir2019{APACrefauthors}Hordoir, R., Axell, L., Höglund, A., Dieterich, C., Fransner, F., Gröger, M.\BDBLHaapala, J. \APACrefYearMonthDay2019. \BBOQ\APACrefatitleNemo-Nordic 1.0: a NEMO-based ocean model for the Baltic and North seas – research and operational applications Nemo-nordic 1.0: a NEMO-based ocean model for the baltic and north seas – research and operational applications.\BBCQ \APACjournalVolNumPagesGeoscientific Model Development121363–386. {APACrefDOI} 10.5194/gmd-12-363-2019 \PrintBackRefs\CurrentBib
Hunter (\APACyear2007) \APACinsertmetastarhunter2007{APACrefauthors}Hunter, J\BPBID. \APACrefYearMonthDay2007. \BBOQ\APACrefatitleMatplotlib: A 2D graphics environment Matplotlib: A 2d graphics environment.\BBCQ \APACjournalVolNumPagesComputing in Science & Engineering9390–95. {APACrefDOI} 10.1109/MCSE.2007.55 \PrintBackRefs\CurrentBib
Huthnance (\APACyear1991) \APACinsertmetastarhuthnance1991{APACrefauthors}Huthnance, J. \APACrefYearMonthDay1991. \BBOQ\APACrefatitlePhysical oceanography of the North Sea Physical oceanography of the north sea.\BBCQ \APACjournalVolNumPagesOcean and Shoreline Management163-4199–231. {APACrefDOI} 10.1016/0951-8312(91)90005-m \PrintBackRefs\CurrentBib
Kalmikov \BBA Heimbach (\APACyear2014) \APACinsertmetastarkalmikov2014{APACrefauthors}Kalmikov, A\BPBIG.\BCBT \BBA Heimbach, P. \APACrefYearMonthDay2014. \BBOQ\APACrefatitleA Hessian-Based Method for Uncertainty Quantification in Global Ocean State Estimation A hessian-based method for uncertainty quantification in global ocean state estimation.\BBCQ \APACjournalVolNumPagesSIAM Journal on Scientific Computing365S267–S295. {APACrefDOI} 10.1137/130925311 \PrintBackRefs\CurrentBib
Kanarik \BOthers. (\APACyear2021) \APACinsertmetastarkanarik2021{APACrefauthors}Kanarik, H., Tuomi, L., Björkqvist, J\BHBIV.\BCBL \BBA Kärnä, T. \APACrefYearMonthDay2021. \BBOQ\APACrefatitleImproving Baltic Sea wave forecasts using modelled surface currents Improving Baltic Sea wave forecasts using modelled surface currents.\BBCQ \APACjournalVolNumPagesOcean Dynamics716-7635–653. {APACrefDOI} 10.1007/s10236-021-01455-y \PrintBackRefs\CurrentBib
Kurapov \BOthers. (\APACyear2011) \APACinsertmetastarkurapov2011{APACrefauthors}Kurapov, A\BPBIL., Foley, D., Strub, P\BPBIT., Egbert, G\BPBID.\BCBL \BBA Allen, J\BPBIS. \APACrefYearMonthDay2011. \BBOQ\APACrefatitleVariational assimilation of satellite observations in a coastal ocean model off Oregon Variational assimilation of satellite observations in a coastal ocean model off Oregon.\BBCQ \APACjournalVolNumPagesJournal of Geophysical Research116C5. \PrintBackRefs\CurrentBib
Kärnä \BOthers. (\APACyear2011) \APACinsertmetastarkarna2011{APACrefauthors}Kärnä, T., de Brye, B., Gourgue, O., Lambrechts, J., Comblen, R., Legat, V.\BCBL \BBA Deleersnijder, E. \APACrefYearMonthDay2011. \BBOQ\APACrefatitleA fully implicit wetting–drying method for DG-FEM shallow water models, with an application to the Scheldt Estuary A fully implicit wetting–drying method for DG-FEM shallow water models, with an application to the Scheldt estuary.\BBCQ \APACjournalVolNumPagesComputer Methods in Applied Mechanics and Engineering2005-8509–524. {APACrefDOI} 10.1016/j.cma.2010.07.001 \PrintBackRefs\CurrentBib
Kärnä \BOthers. (\APACyear2018) \APACinsertmetastarkarna2018{APACrefauthors}Kärnä, T., Kramer, S\BPBIC., Mitchell, L., Ham, D\BPBIA., Piggott, M\BPBID.\BCBL \BBA Baptista, A\BPBIM. \APACrefYearMonthDay2018. \BBOQ\APACrefatitleThetis coastal ocean model: discontinuous Galerkin discretization for the three-dimensional hydrostatic equations Thetis coastal ocean model: discontinuous Galerkin discretization for the three-dimensional hydrostatic equations.\BBCQ \APACjournalVolNumPagesGeoscientific Model Development11114359–4382. {APACrefDOI} 10.5194/gmd-11-4359-2018 \PrintBackRefs\CurrentBib
Kärnä \BOthers. (\APACyear2021) \APACinsertmetastarkarna2021{APACrefauthors}Kärnä, T., Ljungemyr, P., Falahat, S., Ringgaard, I., Axell, L., Korabel, V.\BDBLet al. \APACrefYearMonthDay2021. \BBOQ\APACrefatitleNemo-Nordic 2.0: operational marine forecast model for the Baltic Sea Nemo-nordic 2.0: operational marine forecast model for the Baltic Sea.\BBCQ \APACjournalVolNumPagesGeoscientific Model Development1495731–5749. {APACrefDOI} 10.5194/gmd-14-5731-2021 \PrintBackRefs\CurrentBib
Large \BBA Yeager (\APACyear2008) \APACinsertmetastarlarge2008{APACrefauthors}Large, W\BPBIG.\BCBT \BBA Yeager, S\BPBIG. \APACrefYearMonthDay2008. \BBOQ\APACrefatitleThe global climatology of an interannually varying air–sea flux data set The global climatology of an interannually varying air–sea flux data set.\BBCQ \APACjournalVolNumPagesClimate Dynamics332-3341–364. {APACrefDOI} 10.1007/s00382-008-0441-3 \PrintBackRefs\CurrentBib
Le Dimet \BBA Talagrand (\APACyear1986) \APACinsertmetastarledimet1986{APACrefauthors}Le Dimet, F\BHBIX.\BCBT \BBA Talagrand, O. \APACrefYearMonthDay1986. \BBOQ\APACrefatitleVariational algorithms for analysis and assimilation of meteorological observations: theoretical aspects Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects.\BBCQ \APACjournalVolNumPagesTellus A: Dynamic Meteorology and Oceanography38297–110. \PrintBackRefs\CurrentBib
Leppäranta \BBA Myrberg (\APACyear2009) \APACinsertmetastarlepparanta2009{APACrefauthors}Leppäranta, M.\BCBT \BBA Myrberg, K. \APACrefYear2009. \APACrefbtitlePhysical oceanography of the Baltic Sea Physical oceanography of the baltic sea. \APACaddressPublisherSpringer Science & Business Media. {APACrefDOI} https://doi.org/10.1007/978-3-540-79703-6 \PrintBackRefs\CurrentBib
LeVeque (\APACyear2002) \APACinsertmetastarleveque2002{APACrefauthors}LeVeque, R\BPBIJ. \APACrefYear2002. \APACrefbtitleFinite Volume Methods for Hyperbolic Problems Finite volume methods for hyperbolic problems. \APACaddressPublisherCambridge University Press. {APACrefDOI} 10.1017/CBO9780511791253 \PrintBackRefs\CurrentBib
Loose \BBA Heimbach (\APACyear2021) \APACinsertmetastarloose2021{APACrefauthors}Loose, N.\BCBT \BBA Heimbach, P. \APACrefYearMonthDay2021. \BBOQ\APACrefatitleLeveraging Uncertainty Quantification to Design Ocean Climate Observing Systems Leveraging uncertainty quantification to design ocean climate observing systems.\BBCQ \APACjournalVolNumPagesJournal of Advances in Modeling Earth Systems134. {APACrefDOI} 10.1029/2020ms002386 \PrintBackRefs\CurrentBib
Marotzke \BOthers. (\APACyear1999) \APACinsertmetastarmarotzke1999{APACrefauthors}Marotzke, J., Giering, R., Zhang, K\BPBIQ., Stammer, D., Hill, C.\BCBL \BBA Lee, T. \APACrefYearMonthDay1999. \BBOQ\APACrefatitleConstruction of the adjoint MIT ocean general circulation model and application to Atlantic heat transport sensitivity Construction of the adjoint MIT ocean general circulation model and application to atlantic heat transport sensitivity.\BBCQ \APACjournalVolNumPagesJournal of Geophysical Research: Oceans104C1229529-29547. {APACrefDOI} 10.1029/1999JC900236 \PrintBackRefs\CurrentBib
Maßmann (\APACyear2010) \APACinsertmetastarmassmann2010{APACrefauthors}Maßmann, S. \APACrefYearMonthDay2010. \BBOQ\APACrefatitleSensitivities of an adjoint, unstructured mesh, tidal model on the European Continental Shelf Sensitivities of an adjoint, unstructured mesh, tidal model on the european continental shelf.\BBCQ \APACjournalVolNumPagesOcean Dynamics6061463–1477. {APACrefDOI} 10.1007/s10236-010-0347-6 \PrintBackRefs\CurrentBib
Mellor \BOthers. (\APACyear1998) \APACinsertmetastarmellor1998{APACrefauthors}Mellor, G\BPBIL., Oey, L\BHBIY.\BCBL \BBA Ezer, T. \APACrefYearMonthDay1998. \BBOQ\APACrefatitleSigma Coordinate Pressure Gradient Errors and the Seamount Problem Sigma coordinate pressure gradient errors and the seamount problem.\BBCQ \APACjournalVolNumPagesJournal of Atmospheric and Oceanic Technology1551122–1131. {APACrefDOI} 10.1175/1520-0426(1998)015¡1122:scpgea¿2.0.co;2 \PrintBackRefs\CurrentBib
Met Office (\APACyear2010 - 2022) \APACinsertmetastariris2022{APACrefauthors}Met Office. \APACrefYearMonthDay2010 - 2022. \BBOQ\APACrefatitleIris: A Python package for analysing and visualising meteorological and ocean ographic data sets Iris: A python package for analysing and visualising meteorological and ocean ographic data sets\BBCQ (\PrintOrdinalv3.2 \BEd) [\bibcomputersoftwaremanual]. \APACaddressPublisherExeter, Devon. \PrintBackRefs\CurrentBib
Mohrholz (\APACyear2018) \APACinsertmetastarmohrholz2018{APACrefauthors}Mohrholz, V. \APACrefYearMonthDay2018. \BBOQ\APACrefatitleMajor Baltic Inflow Statistics – Revised Major Baltic inflow statistics – revised.\BBCQ \APACjournalVolNumPagesFrontiers in Marine Science5. {APACrefDOI} 10.3389/fmars.2018.00384 \PrintBackRefs\CurrentBib
Mohrholz \BOthers. (\APACyear2015) \APACinsertmetastarmohrholz2015{APACrefauthors}Mohrholz, V., Naumann, M., Nausch, G., Krüger, S.\BCBL \BBA Gräwe, U. \APACrefYearMonthDay2015. \BBOQ\APACrefatitleFresh oxygen for the Baltic Sea — An exceptional saline inflow after a decade of stagnation Fresh oxygen for the Baltic Sea — an exceptional saline inflow after a decade of stagnation.\BBCQ \APACjournalVolNumPagesJournal of Marine Systems148152–166. {APACrefDOI} 10.1016/j.jmarsys.2015.03.005 \PrintBackRefs\CurrentBib
Moore \BOthers. (\APACyear2004) \APACinsertmetastarmoore2004{APACrefauthors}Moore, A\BPBIM., Arango, H\BPBIG., Lorenzo, E\BPBID., Cornuelle, B\BPBID., Miller, A\BPBIJ.\BCBL \BBA Neilson, D\BPBIJ. \APACrefYearMonthDay2004. \BBOQ\APACrefatitleA comprehensive ocean prediction and analysis system based on the tangent linear and adjoint of a regional ocean model A comprehensive ocean prediction and analysis system based on the tangent linear and adjoint of a regional ocean model.\BBCQ \APACjournalVolNumPagesOcean Modelling71-2227–258. {APACrefDOI} 10.1016/j.ocemod.2003.11.001 \PrintBackRefs\CurrentBib
Rathgeber \BOthers. (\APACyear2016) \APACinsertmetastarrathgeber2016{APACrefauthors}Rathgeber, F., Ham, D\BPBIA., Mitchell, L., Lange, F., Michael a nd Luporini, Mcrae, A\BPBIT\BPBIT., Bercea, G\BHBIT.\BDBLKelly, P\BPBIH\BPBIJ. \APACrefYearMonthDay2016. \BBOQ\APACrefatitleFiredrake: Automating the Finite Element Method by Composing Abstractions Firedrake: Automating the finite element method by composing abstractions.\BBCQ \APACjournalVolNumPagesACM Trans. Math. Softw.43324:1–24:27. {APACrefDOI} 10.1145/2998441 \PrintBackRefs\CurrentBib
Rivière (\APACyear2008) \APACinsertmetastarriviere2008{APACrefauthors}Rivière, B. \APACrefYear2008. \APACrefbtitleDiscontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations: Theory and Implementation Discontinuous Galerkin methods for solving elliptic and parabolic equations: Theory and implementation (\BNUM 35). \APACaddressPublisherSIAM. {APACrefDOI} 10.1137/1.9780898717440 \PrintBackRefs\CurrentBib
Roe (\APACyear1981) \APACinsertmetastarroe1981{APACrefauthors}Roe, P. \APACrefYearMonthDay1981. \BBOQ\APACrefatitleApproximate Riemann solvers, parameter vectors, and difference schemes Approximate riemann solvers, parameter vectors, and difference schemes.\BBCQ \APACjournalVolNumPagesJournal of Computational Physics432357–372. {APACrefDOI} 10.1016/0021-9991(81)90128-5 \PrintBackRefs\CurrentBib
Sirkes \BBA Tziperman (\APACyear1997) \APACinsertmetastarsirkes1997{APACrefauthors}Sirkes, Z.\BCBT \BBA Tziperman, E. \APACrefYearMonthDay1997. \BBOQ\APACrefatitleFinite Difference of Adjoint or Adjoint of Finite Difference? Finite difference of adjoint or adjoint of finite difference?\BBCQ \APACjournalVolNumPagesMonthly Weather Review125123373–3378. {APACrefDOI} 10.1175/1520-0493(1997)125¡3373:fdoaoa¿2.0.co;2 \PrintBackRefs\CurrentBib
Sonnewald \BOthers. (\APACyear2021) \APACinsertmetastarsonnewald2021{APACrefauthors}Sonnewald, M., Lguensat, R., Jones, D\BPBIC., Dueben, P\BPBID., Brajard, J.\BCBL \BBA Balaji, V. \APACrefYearMonthDay2021. \BBOQ\APACrefatitleBridging observations, theory and numerical simulation of the ocean using machine learning Bridging observations, theory and numerical simulation of the ocean using machine learning.\BBCQ \APACjournalVolNumPagesEnvironmental Research Letters167073008. {APACrefDOI} 10.1088/1748-9326/ac0eb0 \PrintBackRefs\CurrentBib
Taylor (\APACyear2001) \APACinsertmetastartaylor2001{APACrefauthors}Taylor, K\BPBIE. \APACrefYearMonthDay2001. \BBOQ\APACrefatitleSummarizing multiple aspects of model performance in a single diagram Summarizing multiple aspects of model performance in a single diagram.\BBCQ \APACjournalVolNumPagesJournal of Geophysical Research: Atmospheres106D77183–7192. {APACrefDOI} 10.1029/2000jd900719 \PrintBackRefs\CurrentBib
Thacker \BBA Long (\APACyear1988) \APACinsertmetastarthacker1988{APACrefauthors}Thacker, W\BPBIC.\BCBT \BBA Long, R\BPBIB. \APACrefYearMonthDay1988. \BBOQ\APACrefatitleFitting dynamics to data Fitting dynamics to data.\BBCQ \APACjournalVolNumPagesJ. Geophys. Res.93C21227. \PrintBackRefs\CurrentBib
Vidard \BOthers. (\APACyear2015) \APACinsertmetastarvidard2015{APACrefauthors}Vidard, A., Bouttier, P\BHBIA.\BCBL \BBA Vigilant, F. \APACrefYearMonthDay2015. \BBOQ\APACrefatitleNEMOTAM: tangent and adjoint models for the ocean modelling platform NEMO NEMOTAM: tangent and adjoint models for the ocean modelling platform NEMO.\BBCQ \APACjournalVolNumPagesGeoscientific Model Development841245–1257. {APACrefDOI} 10.5194/gmd-8-1245-2015 \PrintBackRefs\CurrentBib
Villaret \BOthers. (\APACyear2016) \APACinsertmetastarvillaret2016{APACrefauthors}Villaret, C., Kopmann, R., Wyncoll, D., Riehme, J., Merkel, U.\BCBL \BBA Naumann, U. \APACrefYearMonthDay2016. \BBOQ\APACrefatitleFirst-order uncertainty analysis using Algorithmic Differentiation of morphodynamic models First-order uncertainty analysis using Algorithmic Differentiation of morphodynamic models.\BBCQ \APACjournalVolNumPagesComputers & Geosciences90144–151. \PrintBackRefs\CurrentBib
Virtanen \BOthers. (\APACyear2020) \APACinsertmetastarscipy2020{APACrefauthors}Virtanen, P., Gommers, R., Oliphant, T\BPBIE., Haberland, M., Reddy, T., Cournapeau, D.\BDBLSciPy 1.0 Contributors \APACrefYearMonthDay2020. \BBOQ\APACrefatitleSciPy 1.0: Fundamental Algorithms for Scientific Computing in Python SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.\BBCQ \APACjournalVolNumPagesNature Methods17261–272. {APACrefDOI} 10.1038/s41592-019-0686-2 \PrintBackRefs\CurrentBib
Wallwork \BOthers. (\APACyear2020) \APACinsertmetastarwallwork2020{APACrefauthors}Wallwork, J\BPBIG., Barral, N., Kramer, S\BPBIC., Ham, D\BPBIA.\BCBL \BBA Piggott, M\BPBID. \APACrefYearMonthDay2020. \BBOQ\APACrefatitleGoal-oriented error estimation and mesh adaptation for shallow water modelling Goal-oriented error estimation and mesh adaptation for shallow water modelling.\BBCQ \APACjournalVolNumPagesSpringer Nature Applied Sciences21–11. {APACrefDOI} 10.1007/s42452-020-2745-9 \PrintBackRefs\CurrentBib
Warder \BOthers. (\APACyear2022) \APACinsertmetastarwarder2022{APACrefauthors}Warder, S\BPBIC., Angeloudis, A.\BCBL \BBA Piggott, M\BPBID. \APACrefYearMonthDay2022. \BBOQ\APACrefatitleSedimentological data-driven bottom friction parameter estimation in modelling Bristol Channel tidal dynamics Sedimentological data-driven bottom friction parameter estimation in modelling bristol channel tidal dynamics.\BBCQ \APACjournalVolNumPagesOcean Dynamics. {APACrefDOI} 10.1007/s10236-022-01507-x \PrintBackRefs\CurrentBib
Warder \BOthers. (\APACyear2021) \APACinsertmetastarwarder2021{APACrefauthors}Warder, S\BPBIC., Horsburgh, K\BPBIJ.\BCBL \BBA Piggott, M\BPBID. \APACrefYearMonthDay2021. \BBOQ\APACrefatitleAdjoint-based sensitivity analysis for a numerical storm surge model Adjoint-based sensitivity analysis for a numerical storm surge model.\BBCQ \APACjournalVolNumPagesOcean Modelling160101766. {APACrefDOI} https://doi.org/10.1016/j.ocemod.2021.101766 \PrintBackRefs\CurrentBib
Wunsch \BBA Heimbach (\APACyear2007) \APACinsertmetastarwunsch2007{APACrefauthors}Wunsch, C.\BCBT \BBA Heimbach, P. \APACrefYearMonthDay2007. \BBOQ\APACrefatitlePractical global oceanic state estimation Practical global oceanic state estimation.\BBCQ \APACjournalVolNumPagesPhysica D2301-2197–208. \PrintBackRefs\CurrentBib
Yaremchuk \BOthers. (\APACyear2016) \APACinsertmetastaryaremchuk2016{APACrefauthors}Yaremchuk, M., Martin, P., Koch, A.\BCBL \BBA Beattie, C. \APACrefYearMonthDay2016. \BBOQ\APACrefatitleComparison of the adjoint and adjoint-free 4dVar assimilation of the hydrographic and velocity observations in the Adriatic Sea Comparison of the adjoint and adjoint-free 4dvar assimilation of the hydrographic and velocity observations in the adriatic sea.\BBCQ \APACjournalVolNumPagesOcean Modelling97129–140. \PrintBackRefs\CurrentBib
Zaron \BOthers. (\APACyear2011) \APACinsertmetastarzaron2011{APACrefauthors}Zaron, E\BPBID., Pradal, M\BHBIA., Miller, P\BPBID., Blumberg, A\BPBIF., Georgas, N., Li, W.\BCBL \BBA Cornuelle, J\BPBIM. \APACrefYearMonthDay2011. \BBOQ\APACrefatitleBottom Topography Mapping via Nonlinear Data Assimilation Bottom topography mapping via nonlinear data assimilation.\BBCQ \APACjournalVolNumPagesJournal of Atmospheric and Oceanic Technology28121606–1623. {APACrefDOI} 10.1175/jtech-d-11-00070.1 \PrintBackRefs\CurrentBib
zenodo/Data (\APACyear2022) \APACinsertmetastarzenodo/Data{APACrefauthors}zenodo/Data. \APACrefYearMonthDay2022. \APACrefbtitleThetis Baltic Sea simulation: model and observation data sets [Dataset]. Thetis Baltic Sea simulation: model and observation data sets [Dataset]. {APACrefDOI} 10.5281/zenodo.6497527 \PrintBackRefs\CurrentBib
zenodo/Firedrake (\APACyear2022) \APACinsertmetastarzenodo/Firedrake{APACrefauthors}zenodo/Firedrake. \APACrefYearMonthDay2022. \APACrefbtitleSoftware used in ‘Adjoint-based optimization of a regional water elevation model’ [Software]. Software used in ‘Adjoint-based optimization of a regional water elevation model’ [Software]. {APACrefURL} https://doi.org/10.5281/zenodo.6475963 {APACrefDOI} 10.5281/zenodo.6475963 \PrintBackRefs\CurrentBib
zenodo/Thetis (\APACyear2022) \APACinsertmetastarzenodo/Thetis{APACrefauthors}zenodo/Thetis. \APACrefYearMonthDay2022. \APACrefbtitleThetis coastal ocean model source code [Software]. Thetis coastal ocean model source code [Software]. {APACrefDOI} 10.5281/zenodo.6475560 \PrintBackRefs\CurrentBib
zenodo/Uptide (\APACyear2020) \APACinsertmetastarzenodo/Uptide{APACrefauthors}zenodo/Uptide. \APACrefYearMonthDay2020. \APACrefbtitleFirst release of Uptide v1.0. First release of uptide v1.0. {APACrefDOI} 10.5281/zenodo.3909651 \PrintBackRefs\CurrentBib
Zhang \BOthers. (\APACyear2011) \APACinsertmetastarzhang2011{APACrefauthors}Zhang, J., Lu, X., Wang, P.\BCBL \BBA Wang, Y\BPBIP. \APACrefYearMonthDay2011. \BBOQ\APACrefatitleStudy on linear and nonlinear bottom friction parameterizations for regional tidal models using data assimilation Study on linear and nonlinear bottom friction parameterizations for regional tidal models using data assimilation.\BBCQ \APACjournalVolNumPagesContinental Shelf Research316555–573. {APACrefDOI} 10.1016/j.csr.2010.12.011 \PrintBackRefs\CurrentBib