!--LUA BEGIN ! secondsperyear = 365.25*24.0*3600.0 ! T_ref = 200.0 --- [K] ! CapA = 7.253 --- [J kg-1 K-2] ! CapB = 146.3 --- [J kg-1 K-1] ! ! -- ## Surface temp ----> Surface Enthalpy ! -- See: https://github.com/ElmerCSC/elmerfem/blob/devel/elmerice/Tests/Enthalpy/enthalpy.sif ! -- for reference on boundary conditions ! function surfenthalpy(Y) ! m = -0.01 ! b = 273.15 ! T = m*Y + b ! E = (CapA/2*(T^2 - T_ref^2) + CapB*(T-T_ref)) * secondsperyear^2 ! return E ! end ! ! -- ## Conditional LUA statement ! function IfThenElse(condition,t,f) ! if condition then return t else return f end ! end !---LUA END #yearinsec = 365.25*24*60*60 ![s a^-1] #T_ref = 200.0 ![K] #CapA = 7.253 ![J kg-1 K-2] #CapB = 146.3 ![J kg-1 K-1] Header CHECK KEYWORDS Warn Mesh DB "." "testglacier" Include Path "" Results Directory "" End Constants Stefan Boltzmann = Real 5.67e-08 ![W m−2 K−4] ! Enthalpy Related Params T_ref_enthalpy = real 200.0 ![K] L_heat = real 334000.0 ![J kg-1] ! Cp(T) = A*T + B Enthalpy Heat Capacity A = real 7.253 ![J kg-1 K-2] Enthalpy Heat Capacity B = real 146.3 ![J kg-1 K-1] P_triple = real 0.061173 ![MPa] Triple point pressure for water (Gilbert et al. 2014) P_surf = real 0.1013 ![MPa] Surface atmospheric pressure beta_clapeyron = real 0.0974 ![K MPa-1] clausus clapeyron relationship (Gilbert et al. 2014) End Simulation Max Output Level = 4 Coordinate System = "Cartesian 2D" Simulation Type = "Steady" Steady State Min Iterations = 1 Steady State Max Iterations = 20 Output Intervals = 1 Post File = "diagnostic_enthalpy.vtu" Initialize Dirichlet Conditions = Logical False Restart File = "Stokes_ELA400_diagnostic_temp_constrained.result" Restart Position = 0 End Body 1 Name = "Glacier" Body Force = 1 Equation = 1 Material = 1 Initial Condition = 1 End Equation 1 Name = "Equation1" Convection = "computed" Flow Solution Name = String "Flow Solution" Active Solvers(4) = 1 2 3 4 End Initial Condition 1 ! Pressure = Real 0.0 ! Velocity 1 = Real 0.0 ! Velocity 2 = Real 0.0 ! Velocity 3 = Real 0.0 enthalpy_h = real 125656.0 ![J kg-1] !Temperature = real -3.0 End Solver 1 Exec Solver = "before timestep" Equation = "MapCoordinate_ini" Procedure = "StructuredMeshMapper" "StructuredMeshMapper" Active Coordinate = Integer 2 ! the mesh-update is y-direction Mesh Velocity Variable = String "Mesh Velocity 2" ! The 1st value is special as the mesh velocity could be unrelistically high Mesh Velocity First Zero = Logical True Dot Product Tolerance = Real 0.01 ! Essential to keep a minimum of 10 m flow height Correct Surface = Logical True Minimum Height = Real 10.0 ![m] End Solver 2 Equation = "HeightDepth" !Exec Solver = "Before Timestep" Procedure = "StructuredProjectToPlane" "StructuredProjectToPlane" Active Coordinate = Integer 2 Operator 1 = depth Operator 2 = height End Solver 3 Equation = "Navier-Stokes" Optimize Bandwidth = Logical True Linear System Solver = Direct Linear System Direct Method = "UMFPACK" Linear System Max Iterations = 5000 Linear System Convergence Tolerance = 1.0E-08 Linear System Abort Not Converged = False Linear System Preconditioning = "ILU2" Linear System Residual Output = 10 Steady State Convergence Tolerance = 1.0E-03 ! Stabilization Method = [Stabilized,P2/P1,Bubbles] Stabilization Method = "Bubbles" Nonlinear System Convergence Tolerance = 1.0E-05 Nonlinear System Convergence Measure = Solution Nonlinear System Max Iterations = 25 ! 1 try without non-lin iters Nonlinear System Newton After Iterations = 3 Nonlinear System Newton After Tolerance = 1.0E-03 Nonlinear System Reset Newton = Logical True ! Exported Variable 1 = -dofs 3 "Mesh Velocity" ! Exported Variable 2 = -dofs 3 "Mesh Update" ! Exported Variable 1 = String "Mass Balance" ! Exported Variable 1 DOFs = 1 ! ! Exported Variable 2 = String "Surf Enth" ! Exported Variable 2 DOFs = 1 ! ! Exported Variable 3 = String "Densi" ! Exported Variable 3 DOFs = 1 ! ! Exported Variable 4 = String "Firn" ! Exported Variable 4 DOFs = 1 ! ! Exported Variable 5 = String "Mesh Velocity" ! Exported Variable 5 DOFS = 2 End Solver 4 Equation = String "Enthalpy Equation" Procedure = File "ElmerIceSolvers" "EnthalpySolver" Variable = String "Enthalpy_h" Linear System Solver = "UMFPACK" !Linear System Iterative Method = "BiCGStab" Linear System Max Iterations = 500 Linear System Convergence Tolerance = 1.0E-07 Linear System Abort Not Converged = True Linear System Residual Output = 1 Steady State Convergence Tolerance = 1.0E-04 Nonlinear System Convergence Tolerance = 1.0E-03 Nonlinear System Max Iterations = 10 Nonlinear System Relaxation Factor = Real 1.0 Apply Dirichlet = Logical True Stabilize = True Exported Variable 1 = String "Phase Change Enthalpy" ! (J kg-1) Exported Variable 1 DOFs = 1 Exported Variable 2 = String "Water Content" ! (%) Exported Variable 2 DOFs = 1 Exported Variable 3 = String "temperature" ! (°C) Exported Variable 3 DOFs = 1 End Material 1 Name = "ice" Density = Real #910.0*1.0E-06*(31556926.0)^(-2.0) ! [kg m-3] ---> [MPa a2 m-2] ! Viscosity Model = String "Glen" ! ! Viscosity has to be set to a dummy value to avoid warning output from Elmer ! Viscosity = Real 1.0 ! Glen Exponent = Real 3.0 ! Critical Shear Rate = Real 1.0e-10 ! Rate Factor 1 = Real 1.258e13 ! [MPa^-3a^-1'] ! Rate Factor 2 = Real 6.046e28 ! [MPa^-3a^-1'] ! ! these are in SI units - no problem, as long as ! ! the gas constant also is ! Activation Energy 1 = Real 60.0e3 ! Activation Energy 2 = Real 139.0e3 ! Glen Enhancement Factor = Real 1.0 ! ! the temperature to switch between the two regimes in the flow law ! Limit Temperature = Real -10.0 ! [C] ! Temperature Field Variable = String "temperature" Viscosity Model = String "power law" Viscosity = Real 0.170998e0 ! [MPa a^(1/3)] Viscosity Exponent = Real #1.0/3.0 Critical Shear Rate = Real 1.0e-10 Enthalpy Density = real #910.0 / (1.0e6*yearinsec^2) ! [kg m-3] --> [MPa a2 m-2] Enthalpy Heat Diffusivity = real #(2.1/2050.0) * (1.0e6*yearinsec^-1) ! [kg m-1 s-1] --> [MPa a] Enthalpy Water Diffusivity = real #1.045e-4 / (1.0e6*yearinsec) ! [kg m-1 s-1] --> [MPa a] End !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Body Force 1 Heat Source = real 0.0 Flow BodyForce 1 = Real 0.0 Flow BodyForce 2 = Real #-9.81 * (31556926.0)^(2.0) ! [m a-2] End Boundary Condition 1 Name = "bedrock" Target Boundaries = 1 ! include the bedrock DEM, which has two colums Bottom Surface = Variable Coordinate 1 Real cubic include "steady_ELA400_bedrock.dat" End Normal-Tangential Velocity = True Velocity 1 = Real 0.0e0 Slip Coefficient 2 = Variable Coordinate 2 Real LUA "IfThenElse(((tx[0] > 300) and (tx[0] < 400)), 0.01, 1000.01)" Enthalpy Heat Flux BC = logical True Enthalpy Heat Flux = Real # 0.200 * (31556926.0)*1.0E-06 ! [W m-2] -> [MPa m a-1] End Boundary Condition 2 Name = "sides" Target Boundaries(2) = 3 4 ! combine left and right boundary Velocity 1 = Real 0.0e0 End Boundary Condition 3 Name = "surface" Target Boundaries = 2 ! include the surface DEM, which has two colums Top Surface = Variable Coordinate 1 Real cubic include "steady_ELA400_surface.dat" End Enthalpy_h = Variable Coordinate 2 Real lua "surfenthalpy(tx[0])" !real MATC "3.626*(-2.5-(tx-0.0)*0.006+273.15)^2+146.3*(-2.5-(tx-0.0)*0.006+273.15)-200*(146.3+3.626*200)" !Real lua "-25000.0/150.0*(tx[0]-0)+140000.0" End