{'fraction_non_alphanumeric': 0.03721608348680172, 'fraction_numerical': 0.012622774708410068, 'mean_word_length': 4.569444444444445, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A variational neural network approach for glacier modelling with nonlinear rheology

In this paper, we propose a mesh-free method to solve full stokes equation which models the glacier movement with nonlinear rheology. Our approach is inspired by the Deep-Ritz method proposed in [12]. We first formulate the solution of non-Newtonian ice flow model into the minimizer of a variational integral with boundary constraints. The solution is then approximated by a deep neural network whose loss function is the variational integral plus soft constraint from the mixed boundary conditions. Instead of introducing mesh grids or basis functions to evaluate the loss function, our method only requires uniform samplers of the domain and boundaries. To address instability in real-world scaling, we re-normalize the input of the network at the first layer and balance the regularizing factors for each individual boundary. Finally, we illustrate the performance of our method by several numerical experiments, including a 2D model with analytical solution, Arolla glacier model with real scaling and a 3D model with periodic boundary conditions. Numerical results show that our proposed method is efficient in solving the non-Newtonian mechanics arising from glacier modeling with nonlinear rheology.


Introduction
In recent years, deep neural networks (DNNs) have achieved unprecedented levels of success in a broad range of areas such as computer vision, speech recognition, natural language processing, and health sciences, producing results comparable or superior to human experts typically face a combination of challenges, such as the presence of local features emerging from the nonlinear rheology, saddle point problems due to the incompressibility condition, complex problem domain with high aspect ratios, and high contrast boundary conditions. There has been a lot of effort in developing efficient numerical methods to address these challenges. For example, the work of [5,14,20,23,51,55] designed and analyzed efficient and accurate high-order finite element discretization schemes. Delicated adaptation strategies, see [23,48] and references therein, are developed to ensure that these high-order methods can successfully resolve local features of the flow field and a wide range of length scales. The resulting discretized nonlinear systems have to be solved using either Picard fixedpoint iterations or Newton's method. Most of the recent high-performance solvers adopt modified Newton's iterations together with high-performance iterative linear solvers, see, e.g., [6,15,23,36], to obtain a superlinear convergence rate.
Our work is inspired by the deep Ritz method proposed in [12] and our recent progress in developing deep learning method to solve interface problem [57]. We first formulate the PDEs into a variational problem, which leads to the objective function used by the neural network training. The second part consists of penalty terms arising from the boundary constraint of the governing PDE. In real world models, e.g. glacier sliding, the scale of the problem domain and the scale of physical parameters vary significantly. To address these challenges, we introduce a normalizing layer following the input of neural networks and design a strategy to balance penalty factors due to different boundary conditions. Using this combination of strategies, our proposed deep learning method is capable to solve nonlinear Stokes problems with different geometric scales, different parameter scales and boundary conditions using a universal configuration of network parameters.
The rest of the paper is organized as follows. In Section 2, we introduce the background of the Non-Newtonian ice flow model. In Section 3, we review the basic idea of deep neural network and formulate the variational neural network approach for solving the Non-Newtonian ice flow model. In addition, we also discuss some details of the implementation of our method. In Section 4, we present numerical results to demonstrate the accuracy of our method. Section 5 offers some concluding remarks.

Background
Since glaciers form one of the natural low-pass filters of atmospheric variability, modelling the mechanics of glaciers is instrumental in revealing slow changes in the climate system that might otherwise be obscured by short-term noise. We consider a canonical glacier model that treats the flow of ice as non-Newtonian, viscous, incompressible, and isothermal fluid in the steady-state [11]. In this section, we first present the strong form of the ice flow model, and then discuss its variational formulation that naturally yields the weak form of the ice flow model.

Non-Newtonian ice flow model
For an open, bounded domain Ω ⊂ R d , we denote the velocity field (measured in meters per calendar year, i.e., m a −1 ) and the stress tensor (measured in Pa) of ice flow by u = (u 1 , ..., u d ) : Ω → R d and σ u : Ω → R d×d , respectively. The conservation laws of momentum and mass state that where ρ is the ice density (910 kg m −3 ) and g is gravitational acceleration (9.81 m s −2 , where s represents a second). The stress tensor is split into a deviatoric part and an isotropic pressure p, i.e., σ u = τ u − pI. Denoting the strain rate tensor of the velocity field bẏ the Glen's law of rheology [40] links the deviatoric stress to the strain rate via the constitutive equation where η(u) is the effective viscosity nonlinearly depending on the strain rate and the doubledot product is the Frobenius inner product defined as for two second order tensors σ 1 and σ 2 . Here A = 10 −16 (Pa −n a −1 ) is the flow parameter and we often set n = 3. Using the nonlinear constitutive equation (4), the conservation of momentum (1) becomes −∇ · τ u + ∇p = ρg.
Without loss of generality, we assume the domain Ω is bounded by two disjoint surfaces, Γ b the bottom boundary and Γ t = ∂Ω \ Γ b the top boundary. Let n denote the unit outward normal vector at any point on the boundary ∂Ω. On the top boundary we impose the traction-free boundary condition where p 0 is the atmospheric pressure. Since the atmospheric pressure is often assumed to be negligible [42], so that we adopt p 0 = 0 here. On the bottom boundary we impose a no-penetration condition along the outward normal direction, i.e., and a sliding boundary condition along the tangential direction. Given the map to the tangential direction T = I − n ⊗ n, the sliding boundary condition is given as for some basal drag coefficient β > 0. Note that the basal drag coefficient β is a function of location in general. Since Tn = 0, the sliding boundary condition does not depend on the pressure, which can be reduced to Along the tangential direction, there is a special case that β → ∞, where the sliding boundary condition becomes a non-sliding Dirichlet boundary condition, i.e., Tu = 0. In later numerical examples, our algorithm can be adapted to both cases. To summarize, the ice flow with nonlinear rheology can be modelled by a system of equations where β > 0, f = ρg, T = I − n ⊗ n, and τ u given in (4). With β → ∞, the mixed boundary conditions in (12) and (13) simply become u = 0 on Γ b .

Variational formulation
The momentum equation (9) and the divergence-free mass equation (10) have the variational form: finding a vector function u ∈ H 1 (Ω) satisfying boundary conditions (11) and (12) and a scalar function p ∈ L 2 (Ω) satisfying p = p 0 on Γ t such that for all vector functions v ∈ H 1 (Ω) satisfying v · n = 0 on Γ b and τ v n − p 0 n = 0 on Γ t and all scalar function q ∈ L 2 (Ω) satisfying q = p 0 on Γ t . Applying the divergence theorem for vectors and tensors, the momentum equation (14) can also be expressed as The solutions (u, p) satisfying the above variational form can be expressed as the minimizer of some energy functional subject to the divergence-free constraint (15) and some boundary conditions (11)- (13). This provides a starting point for applying deep neural network approximation to model the glacier flow. We present the energy functional and the corresponding optimization problem for the sliding and the the non-sliding boundary conditions as follows.
Sliding bottom boundary. Subject to the traction free boundary condition (11) on the top boundary, the boundary integral at the top boundary is eliminated, i.e., where p 0 is the atmospheric pressure at the top boundary. At the bottom boundary, applying the condition v · n = 0 on Γ b and the sliding boundary condition (13) we have as range(T) and range(I − T) are mutually orthogonal and (I − T)v = 0. Therefore, the boundary integral in (16) can be simplified as As shown in [11], the solutions (u, p) satisfying the variational form (16) can be considered as the minimizer of the the energy functional subject to the divergence-free condition (10). Formulating this constrained optimization problem using the method of Lagrange multiplier, we have the Lagrangian functional in which the pressure function p plays the role of the Lagrange multiplier. Since the direc- The solution to the Lagrange multiplier, (u, p) such that L s (u, p)[v, q] = 0 for all (v, q) is equivalent to the solution to the variational form defined in (14) and (15). Thus, the velocity field u of the ice flow model can be obtained by minimizing the energy functional E s (u) subject to the divergence-free constraint, the traction-free boundary condition (11), and the Dirichlet boundary conditions (12) and (13).
Non-sliding bottom boundary. With non-sliding boundary condition at the bottom, i.e., u = 0 on Γ b , we do not need to impose the boundary condition (13) using the boundary integral as the sliding boundary case in (19). Following a similar derivation as above, the velocity field u of the ice flow model with non-sliding boundary can be obtained by minimizing the energy functional subject to the divergence-free constraint, the traction-free boundary condition (11), and the Dirichlet boundary conditions u = 0 on Γ b . The corresponding Lagrangian functional becomes

Divergence-free solutions spaces
Instead of searching for the saddle point solution of the Lagrangian functionals, we can construct divergence-free solution spaces for the velocity field u to reduce the constrained optimization problems to unconstrained optimization problems. In the two dimensional case, i.e., (x, y) ∈ Ω ⊂ R 2 , we can define a potential function φ : Ω → R, which leads to a vector function that satisfies the divergence-free condition by construction. In the three dimensional case, i.e., (x, y, z) ∈ Ω ⊂ R 3 , we can define a vector function φ : Ω → R 3 , so that satisfies the divergence-free condition. Using the divergence-free construction of the velocity field, we can directly minimize the energy functional E s (u) and E ns (u) subject to appropriate boundary conditions without using the Lagrangian formulation.

Formulation of the DNN method
We first discuss the background of DNNs and then develop the DNN-based methods for solving the Non-Newtonian ice flow problems.

DNNs and its approximation property
There are two ingredients in defining a DNN. The first one is a linear map of the form The second one is a nonlinear activation function σ : R → R. Common examples of the activation function include the rectified linear unit (ReLU), σ(x) = max(0, x), the soft-plus function, σ(x) = log(e x + 1), the sigmoid function, σ(x) = (1 + e −x ) −1 , etc. The definition of activation function can be trivially extended to a nonlinear map σ : R n → R n by applying the scalar-valued activation function element-wise to each of its inputs.
Using above definitions, we are able to define a continuous function F (x) by a composition of linear transforms and activation maps, i.e., where  (23) as θ ∈ Θ, where θ is a high-dimensional vector and Θ is the space of θ. The DNN representation of a continuous function can be viewed as We use F = {F (·, θ)|θ ∈ Θ} to denote the set of all expressible functions by the DNN parameterized by θ ∈ Θ. Then F provides an efficient way to represent unknown continuous functions, comparing with a linear solution space used in classic numerical methods, e.g., a trial space spaced by linear nodal basis functions in the FEM. In the sequel, we shall discuss the approximation property of the DNN, which is relevant to the study of the expressive power of a DNN model [9,52]. Early investigations of the approximation property of neural networks can be found in [10,21], where the authors analyzed the approximation property for the function class given by a feed-forward neural network with a single hidden layer. Later, many authors analyzed the error estimates of such neural networks in terms of the number of neurons, the number of network layers, and the choice of activation function; see [13,45] for a throughout review of relevant works.
In recent years, DNNs have shown successful applications in a broad range of problems, including classification of complex systems and construction of response surfaces for highdimensional models. Significant efforts have been devoted to investigate the dependence of the expressive power of DNNs on the network configuration. For example, in [9], the authors proved that convolutional DNNs were able to express multivariate functions given in so-called hierarchical tensor formats. In [60], the author studied the expressive power of shallow and deep neural networks with piece-wise linear activation functions and established new rigorous upper and lower bounds for the network complexity in approximating Sobolev spaces. In [18], the authors proved that the continuous piece-wise bilinear finite element space is embedded in the space of DNNs constructed from the ReLU activation function, a sufficient width and a sufficient depth. Thus, we can uses DNNs to approximate the solution space spanned by the FEM basis.

Formulation of the variational neural network approach
The model problem of the ice flow with nonlinear rheology (13) can be solved by using numerical methods such as FEMs and FDMs. However, the nonlinearity in the problem bring essential difficulty in solving (13) by the traditional numerical methods. Inspired by the recent development of deep learning based numerical method for solving variational problems [12,57], we will develop a variational DNN method to solve the glacier model problem defined in (13).
According to the divergence free property, i.e. ∇ · u = 0, we can represent the velocity field u as follows: where φ i is a scalar-valued function approximated by F i (x, θ). Here, F i (x, θ) is a DNN representation with x ∈ R d as its input and scalar output defined in (23). Moreover, θ denotes all the parameters that will be determined during the training stage. Denoting the DNN representation of the velocity field byũ, the numerical solution of (13) can be obtained by finding θ ∈ Θ that minimizes the energy functional E s (ũ) subject to various boundary conditions. Ideally, we want to defineũ using θ ∈ Θ b ⊂ Θ, where Θ b is the maximal parameter subset such that the resultingũ satisfies the traction-free boundary condition (11), and the Dirichlet boundary conditions (12) and sliding bottom boundary condition (13). After parameterizing the expressible function space by θ ∈ Θ b , we equivalently define the variational problem (19) as (26) Note that the divergence-free condition in (15) is automatically satisfied by the representation (25), and thus no additional treatment is needed.
The variational problems (26) is not convex in general and the integrals in (26) do not have a closed-form expression. Thus we numerically approximate the integrals by the Monte Carlo method and use the stochastic gradient descent (SGD) method [4] to minimize the objective function after Monte Carlo approximation. Denoting θ k the kth component of the high-dimensional vector θ, the derivative of J(θ) with respect to θ k can be approximated as After approximating the gradient of J(θ), we can update each component of θ as follows: where η is the learning rate. To accelerate the training of the neural network, we use the Adam version of the SGD method [27].
In the objective function (26), it is rather challenging to restrict the neural network parameter θ to the subset Θ b that satisfies the boundary conditions, because the boundary of the subset may have complicated geometry. To address this issue, we adopt a relaxation approach by imposing boundary conditions as penalty terms. For each of the boundary conditions in (11)-(13), we define a linear map B j that maps the solutionũ to the residual of the boundary constraint, which defines a boundary integral for the boundary ∂D j . Then, the boundary conditions can be imposed as soft constraints via penalty terms to the objective functional J(·) in (26). This leads to a new objective functionalũ where b is the number of boundary conditions needs to be imposed. Note that when a penalty term ε −1 j B j (θ) approaches zero as ε j → 0, the corresponding boundary condition is satisfied in the weak sense.

Implementation details
In addition to the conventional DNN formulation for solving PDE, we also made some modification that are shown to be essential in resolving the scale of the problem domain and the scale of physical parameters in real world problems.
Normalizing Layer. In some glacier model problems, periodic boundary conditions are used. For example, periodic boundary conditions are imposed in the two horizontal directions in the in the 2D model presented in Example C of the benchmarks [42]. To this end, we add an extra layer between input and first dense layer of our network. That maps (x, y) to terms like cos(2πx/T ), sin(2πx/T ), cos(2πy/T ), sin(2πy/T ) , where T is the period in the horizontal dimensions.
In most of the glacier modelling problems, the vertical scaling can be two or three orders of magnitude smaller than the horizontal scaling. To resolve this scaling, we also introduce a reparametrization: x → x/L, where L is the diameter of the domain in corresponding dimension, to the normalization layer from input. Resolving this scaling is essential to ensure the expressibility of neural networks for the problems we are targeting here, because nonlinearity can be represented by a neural network may only exist in a bounded input domain. Taking a scalar single layer perception with the Sigmoid activation function σ(z) = (1 + e −z ) −1 as an example, we have becomes a constant approximately. Thus, to sufficiently express nonlinear functions, scale of W 1 has to be O(L −1 ) without any reparametrization. However, most of existing implementations of machine learning toolboxes only support single-precision float-point numbers. This leads to inaccurate representation of W 1 and the gradient of the objective function, as the length scale of most of the real-world problems can be O(10 4 ) − O(10 6 ). The reparametrization is critical to ensure the expressibility of neural networks and the numerical stability in these real-world problems.
Balancing penalty parameters of boundary conditions. Except for periodic boundary conditions, as shown in (29), the other boundary conditions are enforced by penalty terms. However, there is no general framework to choose the penalty weight ε j in the context of DNN. It may not be suitable to use the same penalty parameter ε j for different boundary conditions, as the magnitude of the boundary maps B j (θ) can vary significantly. For example, in the 3D model presented in Section 4.3, the term B j (θ) is O(10 −4 ) for the traction-free condition at the top boundary and O(10 6 ) for the basal friction condition at the bottom boundary. Thus, finding the correct penalty parameter ε j for each of the boundary conditions is critical in balancing the soft constraints due to various boundary conditions. We propose to first estimate the scaling relations between each of the boundary integrals B j (·) with the objective functional J(·) at the initialization stage of the network training. Denoting the initial parameter by θ 0 , we can then assign ε j as so that the objective functional J(·) and each of the boundary integrals B j (·) can be approximately balanced during the training. Here the common factor ε 0 is chosen to be 50, which is an empirical constant shown to be sufficient in previous research [24] and in our numerical experiments on simpler nonlinear models.
Network Hyperparameters. In the numerical experiments of this work, we will apply the same network to represent φ in both two and three dimensions. After the normalizing layer, it has six latent dense layers with width 10. Between these dense layers we apply the Sigmoid activation function to guarantee the smoothness of our representation. There is no activation function between the second last layer and the output layer. We illustrate the structure of our network in Fig. 1.

Numerical Results
In this section, we shall present numerical results in solving non-Newtonian Stokes equations to demonstrate the performance of the propose method. First we consider a two dimensional synthetic model with an analytical solution to demonstrate the efficiency and accuracy of our method. Then we apply the developed method on realistic benchmark prob-lems [42] to demonstrate its generality. This includes a two dimensional model of the Arolla glacier and a three dimensional box model.

A 2D Model on irregular domain
We start with a two dimensional model defined on a domain that is enclosed by To setup the benchmark, we start with a ground truth potential function which leads to an analytical expression of the velocity field The above velocity field satisfies the conservation of mass condition in Eq. (10) and the no-penetration boundary condition in Eq. (12) at the bottom boundary Γ b . Then, by substituting this velocity field into Eq. (9), we obtain an explicit expression of the forcing term f in Eq. (9), a function p 0 for the traction-free boundary condition Eq. (11) at the top boundary Γ t , and a function β for the sliding boundary condition Eq. (13) at the bottom boundary Γ b .
Since the shape of the domain is a semicircle in this example, it is not necessary to use normalization layers or balance regularization factors. We use a constant regularization factor 1 ε = 50 for all boundary conditions. During training, in every step of stochastic gradient descent, we use 5000 random samples uniformly distributed in the domain Ω to evaluate Eq. (27) and 1000 uniform random samples on the boundary Γ to construct the soft constraints in Eq. (11)-Eq. (13). The learning rate, which is the step length of SGD, is set to be 10 −3 . We renew the training data every 200 steps of learning. The configuration of the neural network is illustrated in Fig. 1.
In Fig. 2, we compare the results of our proposed method with the ground truth. Since the velocity field is invariant to the potential function φ up to a constant shift, we shift the output of φ learned by neural networks to ensure it has the same average as the ground truth in the comparison. For both of the velocity field and the potential function, the neural network offers comparable results with the ground truth.
In Fig. 3 we compare the relative L 2 error of the velocity field with the energy functional E s (u) in Eq. (19). We observe that E s (u) follows a similar trend to the relative L 2 error. The relative L 2 error drops to 0.005 after 10000 steps of training.

The Arolla model
The two-dimensional model of the Arolla glacier is a diagnostic experiment along the central flow-line of a temperate glacier in the European Alps (Haut Glacier d'Arolla, Switzerland), based on earlier experiments by [3,41]. The domain of this model is enclosed by the longitudinal surface and bedrock profiles of Haut Glacier d'Arolla (as shown in Fig. 4 (b) and (c)). This model has been used as a benchmark example in various investigations of the forward modelling and the inverse modelling of glaciers, see [43]   where θ = π/1800. In this example, periodic boundary conditions are applied to boundaries in the horizontal directions, the traction-free boundary conditions in Eq. (11) is applied to the top boundary, and the sliding boundary condition in Eq. (12) and Eq. (13) is prescribed at the bottom boundary. The basal sliding coefficient is defined as β(x, y) = 1000 + 1000 sin 2πx L sin 2πy L .
In Fig.6(b) we plot the value of basal friction coefficient β. Empirically speaking, the surface speed becomes larger when β is small. The ground truth solution is unknown. In Fig.7 we plot the surface velocity against y direction at the slice x = L 4 . In addition we compare our prediction with various result in [42] (experiment C, L = 20km).

Conclusions
In this paper, we investigate deep-learning methods to simulate non-Newtonian ice flow models. By formulating the model problem into a variational form and constructing a divergence-free solution space using DNN, we convert the solution of the nonlinear PDE into the solution of an optimization problem. Despite the rather high-dimensionality of the parameter space of the DNN, the solution can still be efficiently obtained by the SGD method.
In this framework, once we are able to sample from the problem domain and its boundary, no special treatment is needed to handle irregular boundaries of the problem domain. Therefore, the proposed method is easy to implement and mesh-free. To address the application in realworld computation, we introduce normalizing layers and adaptive boundary penalties in the configuration of our network. Finally, we present numerical experiments to demonstrate the performance of the proposed method. Specifically, we use the DNN method to solve ice flow model for the 2D Arrola Glacier and a 3D box model with real-world scaling. Our numerical experiments demonstrate that the DNN method provides satisfactory simulation results for glacier flows. We expect this method can be applied to a more general class of problems in non-Newtonian mechanics that satisfies certain variational principles.