(2 / 3) ϵ pl , dev ij ϵ pl , dev ij . The changes in strength which result from the competing processes of shear-induced bond failure and strength recovery due to rapid sintering processes are described by the phenomenological equations In these equations, Σ$_{R}$ is the residual strength associated with intergranular friction which remains even when all bonds between ice granules have failed, Σ$_{S}$ is the structural strength contribution of intergranular bonds, ε$_{S}$ is a characteristic shear strain associated with bond failure, Σ 0 S is the limit strength reached after prolonged ageing, and τ is the characteristic time constant for strength recovery by sintering/ageing processes. Owing to the possibility of structural softening, deformation may proceed in a localized manner even if the external loading is homogeneous, and spontaneous strain localization is indeed manifest in our experimental data. In such situations, the mathematical formulation of the deformation problem may become ill-posed, leading to mesh dependence of numerical solutions. To mitigate this problem, we adopt the suggestion of Aifantis [15] of adding a second-order gradient of the plastic shear strain to the yield function, a method which has been proven to be thermodynamically consistent by Gurtin and Anand [16] (for alternative regularization methods, see [17]). This leads to the constitutive equation Mathematically, the non-local term suppresses strain localization on scales below the characteristic length l and defines a finite, mesh independent width of the deformation bands. The Laplacian in the yield condition necessitates a corresponding higher-order boundary condition, which we take to be n . ∇ ε pl S = 0 , i.e., we require the gradient of the shear strain to be zero in the direction n perpendicular to the specimen surface. This choice of the higher-order boundary condition is found to correctly reproduce the behavior observed when a deformation band reaches the end of the specimen. In applying our model to the compression experiments considered here, we assume the parameters ν = 0 , α = 1 / 2 . This allows us to reduce the model to three coupled equations for the axial compressive stress s := σ$_{xx}$ where x is the coordinate along the compression axis, the axial strain e := ε pl xx < 0 , and for the structural strength Σ$_{S}$ . A derivation of the simplified equations and a demonstration of the correctness of the assumptions ν ≈ 0 , α ≈ 1 / 2 is given in the Supplementary Material. The axial stress is found to be s = E$_{eff}$ ( e$_{ext}$ − 〈 e 〉 ) , with E − 1 eff = 〈 E − $^{1}$( e ( x )) 〉 (7) S = (Σ$_{R}$ + Σ$_{S}$) exp[ − aε pl V ] + E$_{0}$l $^{2}$Δ ε pl S . (6) Σ = Σ$_{R}$ + Σ$_{S}$ , ∂$_{t}$ Σ$_{S}$ = − ˙ ε pl S ε$_{S}$ Σ$_{S}$ + 1 τ (Σ 0 S − Σ$_{S}$) . (5) 5