Evaluating Simplifications of Subsurface Process 1 Representations for Field-scale Permafrost Hydrology Models

. Permafrost degradation within a warming climate poses a significant environmental 6 threat through both the permafrost carbon feedback and damage to human communities and 7 infrastructure. Understanding this threat relies on better understanding and numerical 8 representation of thermo-hydrological permafrost processes, and the subsequent accurate 9 prediction of permafrost dynamics. All models include simplified assumptions, implying a tradeoff 10 between model complexity and prediction accuracy. The main purpose of this work is to 11 investigate this tradeoff when applying the following commonly made assumptions: (1) assuming 12 equal density of ice and liquid water in frozen soil; (2) neglecting the effect of cryosuction in 13 unsaturated freezing soil; and (3) neglecting advective heat transport during soil freezing and thaw. 14 This study designed a set of 62 numerical experiments using the Advanced Terrestrial Simulator 15 (ATS v1.2) to evaluate the effects of these choices on permafrost hydrological outputs, including 16 both integrated and pointwise quantities. Simulations were conducted under different climate 17 conditions and soil properties from three different sites in both column- and hillslope-scale 18 configurations. Results showed that amongst the three physical assumptions, soil cryosuction is 19 the most crucial yet commonly ignored process. Neglecting cryosuction, on average, can cause 10% 20 ~ 20% error in predicting evaporation, 50% ~ 60% error in discharge, 10% ~ 30% error in thaw 21 depth, and 10% ~ 30% error in soil temperature at 1 m beneath surface. The prediction error for 22 subsurface temperature and water saturation is more obvious at hillslope scales due to the presence 23 of lateral flux. By comparison, using equal ice-liquid density has a minor impact on most 24 hydrological variables, but significantly affects soil water saturation with an averaged 5% ~ 15% 25 error. Neglecting advective heat transport presents the least error, 5% or even much lower, in most 26 variables for a general Arctic tundra system, and can decrease the simulation time at hillslope scales by 40% ~ 80%. By challenging these commonly made assumptions, this work provides permafrost hydrology modelers important context for better choosing the appropriate process

gas to the atmosphere and exacerbating global warming. In Arctic tundra, permafrost also plays 48 an important role in maintaining water, habitat of wildlife, landscape, and infrastructure (Berteaux 49 Walvoord and Kurylyk, 2016). Therefore, the occurrent 53 and potential impacts motivate the development of computational models with the goal of better 54 understanding the thermal and hydrological processes in permafrost regions, and consequently 55 predict permafrost thaw accurately. 56 Simulating soil freezing and thaw processes is a challenging task that incorporates mass and energy 57 transfer among atmosphere, snowpack, land surface (perhaps with free water), and a variably 58 Even in the most process-rich models of permafrost change, three such physics simplifications are 73 often made: representing water at constant density (thereby neglecting the expansion of ice relative 74 to liquid water), neglecting cryosuction of water in unsaturated, partially frozen soils, and 75 neglecting advective heat transport. 76 First, because of the lower density of ice than liquid water, freezing water must expand the volume 77 of the porous media, push liquid water into nearby volume, or otherwise expand the volume 78 occupied by that water. As all of the current set of models operate under the assumption of a rigid 79 solid matrix and thus the absence of mechanical equations describing matrix deformation or frost 80 heave, including this expansion typically results in large pressures that must be offset by grain 81 compressibility or another mechanism. Therefore, the densities of ice and liquid water are Arcos multiphysics library leveraged in ATS; this allows the precise model physics to be specified 167 and configured at runtime through the use of a dependency graph describing swappable 168 components in the model physics . 169

Ice density 170
The density of ice (kg/m 3 ) is represented as a Taylor series expansion in both temperature and 171 pressure: 172 173 and the density of liquid water (kg/m 3 ) is represented as: 174

Cryosuction
where sn is the saturation of n-phase and the subscripts n = l, i, g are liquid, ice, and gas phases, 184 respectively; is a coefficient; Lf is the heat fusion of ice; pn (n = l, g) is the pressure of n-phase; 185 S * is the Van Genuchten model. This physically derived formulation can describe the change of 186 matric suction in the frozen zone due to the change of ice content, and thus has the capacity to 187 represent cryosuction. 188 Alternatively, to exclude the effect of cryosuction in this study, we used the Van Genuchten model 189 to determine the total water content, including liquid water and ice. The liquid water content is In hillslope scenarios, hillslopes with northern and southern aspects are considered to investigate 226 physics representation comparisons under the same climate and soil condition (i.e., at a given site) 227 but different solar radiation incidence. Furthermore, hillslopes with both convergent and divergent 228 geometries are included to compare the sensitivity of simulated discharge on process 229 representation. These scenarios can incorporate many types of Arctic systems at the described plot-230 to-regional scales, but explicitly ignore the effects of microtopography or other local-scale 231 focusing mechanisms such as water tracts or thermo-erosion gullies. The objective is to reach a 232 conclusion on the influence of the three physics representation that can be widely applicable in 233 many Arctic systems. can impede the heat transport between the air and the underlying mineral soil, resulting in varying 279 thaw depth (or permafrost table depth) along a hillslope, which has been observed at the site Teller 280 (Jafarov et al., 2018). In this paper, hillslope meshes were constructed following this observation 281 so that the organic layers are thicker at the top and bottom of a hillslope, as described in the next 282 section.  Table 2. 296

298
In the hillslope scenario, we designed the mesh based on observations at Teller to represent a 299 generalized, varying-thickness low Arctic hillslope. A hillslope mesh was created first by 300 generating a pseudo-2D surface mesh with 50 cells and then extruding the 2D mesh downward by 301 50 m. The pseudo-2D surface was designed in a trapezoidal shape with a single, variable-width 302 cell in the cross-slope direction to represent convergent/divergent hillslopes, the short and long 303 sides of which are 200 m and 800 m, respectively (see Figure 3). Vertically, from surface 304 downward, the grid size distribution was the same as the column mesh for each site. The domain 305 is also composed of three layers, same as the column, while the numbers of cells representing each 306 soil layer (i.e., soil layer thickness) are different along the hillslope. The thickness distribution of 307 the first two layers of each site is shown in Table 3. The third layer of a hillslope for all sites is the 308 principal mineral soil. Additionally, hillslope meshes with different aspects (i.e., north, south) were 309 also created.

Model setup 315
To study how the representations of the three physical processes (i.e., ice density, cryosuction, and 316 advective heat transport) affect simulated hydrological outputs at different scales and hillslope 317 topography features, and under various forcing and soil conditions, 62 model simulations were 318 conducted, summarized in Table 4. To examine the validity of the assumption of equal density 319 between ice and liquid, we included cryosuction and advective heat transport in models. To 320 investigate the role of cryosuction in permafrost modeling, we used different density, while 321 neglecting advective heat transport to decrease the computation cost. Note that neglecting 322 advective heat transport in these runs can reduce the effect of cryosuction on simulation predictions, 323 as cryosuction moves water which would itself advect energy. To compare the difference between 324 neglecting and including heat advection, we used different density expressions for ice and liquid, 325 and included cryosuction. Particularly, in order to understand the impact of advective heat 326 transport on permafrost process when soil is at its wettest, we designed two extreme cases under 327 the warm, wet conditions of the Teller site in which soil evaporation was artificially reduced. 328 These runs were designed to maximize water flux and therefore maximize the potential for 329 advective heat transport to affect predictions. 330

332
Prior to simulating all cases, two steps of initialization are carried out for each site. First, a column 333 model with a given initial water table depth and above-0 ºC temperature was frozen by setting the 334 bottom temperature at a constant value of -10 ºC until a steady-state frozen soil column is formed. 335 The initial water table depth is chosen to ensure that the frozen column's water table, after 336 accounting for expansion of ice, is just below the soil surface. The pressure and temperature 337 profiles of the frozen column were used as the initial conditions of the second step initialization. 338 Before proceeding, the observed forcing data (period of 2011-2020) was averaged across the years 339 to form a one-year, "typical" forcing year, which was then repeated 10 times. Using this typical 340 forcing data and the solutions of the first step, we solved the column model in a transient solution, 341 calculating an annual cyclic steady-state and obtaining the pressure and temperature fields at the 342 end of the 10 th year. The final state was then used as initial condition in formal simulations listed 343 in Table 4. The temperature at bottom was constant at -10 ºC. 344

Evaluation metrics 345
To fully assess the effect of representation of ice density, advective heat transport, and cryosuction 346 in permafrost hydrology modeling, we used root mean squared error (RMSE) and normalized 347 Nash-Sutcliffe efficiency (NNSE) as performance metrics. RMSE has the same dimension with 348 the corresponding variables, which can be used to evaluate the average absolute deviation from a 349 benchmark, defined by: 350 where xt and yt are the two modeled datasets to compare from the initial time point (t = 1) to the 352 end (t = N). 353 NNSE is a normalized dimensionless metric describing the relative relationship between an 354 estimation and a reference, which is oftentimes used for evaluating hydrological models. 355 where the modeled results xt (obtained without physics simplification) is considered as the 357 benchmark, and ̅ is the mean value of the benchmark. NNSE approaching to 1 indicates perfect 358 correspondence between two observations. 359 In addition, we also used normalized mean absolute error (MAE) to quantify the percentage change For variables with zero as the smallest value, such as evaporation, discharge, and thaw depth, the 368 corresponding average value was selected as the reference. 369

Results 370
This section compares simulated outputs over the period of 2011-2020 for the three physics under 371 different simulating conditions. We focus on the impact on integrated variables, such as 372 evaporation, discharge, averaged thaw depth, and depth-dependent variables, such as temperature, 373 and total water saturation (ice and liquid). For hillslope models, we chose five surface locations 374 according to the slope geometry to collect simulated data, which were then averaged to obtain a 375 single outcome for each variable of interest. 376

Ice density 377
To evaluate the representation of ice density on permafrost process simulation, we compared 378 evaporation, discharge, thaw depth, and total water saturation between simulations using equal and 379 different ice density expressions.

408
In addition, we investigated how much the assumption of equal ice-liquid density can affect 409 simulation time at hillslope scale. Using 10-year simulations with real ice density as references, 410 the percentage change of time consumed after applying equal ice-liquid density was calculated and 411 displayed in Figure 6. Overall, under the density assumption, it may take less time (positive 412 percentage), but no more than 25% and on average lower than 10%. However, it may also increase

Cryosuction 421
To evaluate the effect of cryosuction on permafrost process predictions, we compared evaporation, 422 discharge, thaw depth, total water saturation, and temperature obtained through simulations 423 including and neglecting cryosuction.  434 Figure 8 shows the effect of cryosuction on column-scale simulated thaw depth and total water 435 saturation at 5 cm beneath surface. RMSEs of thaw depth range from 3 cm to 8 cm. Though still 436 one order of magnitude smaller than the average annual thaw depth, the estimation error due to 437 neglecting cryosuction is obvious in summer, especially at areas with cold temperature and low 438 rainfall like Barrow. By comparison, at Teller, where the largest thaw depth is over double of 439 Barrow and Sag due to its higher temperature, soil cryosuction does not essentially affect thaw 440 depth compared to the other two sites. Similarly, for the total water saturation, at Barrow, the effect 441 of cryosuction is more clearly observed, not only during cold seasons as observed for density 442 representation (section 4.1), but also in summers. 443

462
Neglecting soil cryosuction has a similar impact on hydrological outputs in hillslope scale models. 463 Figure 10 shows the comparison of the variables discussed above under the Sag climate. 464 Evaporation, thaw depth, and temperature are presented based on south-facing divergent hillslope 465 models, while discharge and water saturation are from hillslope models with north-facing 466 convergent geometry. In general, neglecting soil cryosuction has a smaller effect on integrated 467 variables (evaporation and discharge) compared with other pointwise variables. Though thaw 468 depth presents a high NNSE, approximately 0.94, and low RMSE, about 4.5 cm compared to the 469 average, indicating a good match between models considered and excluded cryosuction, the 470 estimation error during summer may reach as high as 10 cm, particularly from 2011 to 2017, as 471 shown in Figure 10 (c). Obvious errors in water saturation and temperature, similar with column-472 scale models, occur almost annually with respect to extrema during winter and summer. Overall, 473 compared to column-scale models, differences in evaporation, discharge, thaw depth, and surface 474 temperature due to neglecting cryosuction effect are relatively reduced at hillslope scale if 475 comparing NNSEs (Table 5). Localized subsurface variables, such as water saturation and 1m-476 depth soil temperature, show increased errors from column to hillslope scale models, which is 477 primarily caused by lateral flux exchange captured by hillslope modeling. 478

Advective heat transport 486
This section evaluates the performance of advective heat transport in modeling permafrost process. 487 As above, we investigated the influence of neglecting heat advection on evaporation, discharge, 488 thaw depth, total water saturation, and temperature. Overall, heat advection does not play a vital 489 role in a normal Arctic system after comparing all hydrological outputs from models with different at Teller. This specific scenario is chosen to maximize the potential effect of advective heat 497 transport in a hillslope-scale Arctic system. Figure 11 illustrates comparisons on all outputs 498 mentioned above from hillslope models without heat advection and with full thermal 499 representation. Apparently, all RMSEs are extremely small, at least two orders of magnitude lower 500 than the corresponding variable average. Almost all NNSEs are approximately one, even for thaw 501 depth, localized water saturation, and temperature. Therefore, for most Arctic systems at this scale, 502 it is reasonable to neglect advective heat transport. 503

507
In addition to simulated results, we also compared simulation times in percentage change between 508 hillslope models neglecting and including heat advection. ATS uses Algebraic Multigrid method 509 as preconditioner for solving, which has a relatively deficient performance in dealing with 510 hyperbolic equations. Thus, incorporating advective heat transport will aggravate computational 511 cost, particularly in case of both large spatial and temporal scale. Figure 12

Comprehensive comparison 521
In the above three sections, we discussed time-series simulation comparisons. This section will 522 analyze the effect of equal ice-liquid density, neglecting cryosuction, and neglecting heat 523 advection on permafrost modeling outputs from holistic, average perspectives. 524 First, we extracted NNSEs of all variables obtained from all comparing models for qualitative 525 analysis. Table 6 shows an example based on column-scale models under conditions of three 526 different sites. Red numbers highlight the obviously reduced NNSEs of one or two processes 527 among the three for each variable. Overall, neglecting advective heat transport has the least 528 influence on model outputs. Equal ice-liquid density primarily affects saturation and has less effect 529 on other variables. Excluding soil cryosuction makes the greatest impact on almost all variables, 530 especially in a relatively wet environment. Among the variables, evaporation and surface 531 temperature are less affected by the three physical process representations, while location-based 532 water saturation is most affected. 533 * sw and T in Table 6 are water saturation and temperature, respectively.

535
Furthermore, to compare across the physics quantitively, we calculated the mean absolute error 536 (MAE) for each variable of interest over the simulation period of 2011-2020. For evaporation, 537 discharge, and thaw depth, the MAEs are normalized by the corresponding variable average 538 (numbers in Figure 13 (a)); for water saturation and temperature, the MAEs are normalized by 539 their average annual fluctuation range (numbers in Figure 13  From the perspective of 10-year average, in general, each physical process of Arctic system 546 discussed in this paper presents a similar impact on variables between column and hillslope scales. 547 Under climate and soil conditions of three different sites, neglecting cryosuction in permafrost 548 modeling leads to the greatest influence on hydrological prediction amongst the three physical 549 assumptions. As seen in Figure 13 (a), it will result in 10% ~ 20% deviation in evaporation, 50% 550 ~ 60% in discharge, and 10% ~ 30% in thaw depth. Evaporation is the least affected among the 551 three variables. Discharge is more affected in regions with abundant rainfall (Teller), while in 552 regions with less precipitation, evaporation and thaw depth are relatively affected (Barrow). By 553 comparison, assuming equal ice-liquid density and neglecting advective heat transport may only 554 cause 10% and 5% or even much lower error, respectively, in reference to the annual average of a 555 variable. Specially in Barrow, models utilizing the same ice and liquid densities and ignoring 556 advective heat seem to make an obvious impact on discharge, whereas this also results from its 557 extremely low discharge (Figure 7 (b)). 558 Figure 13 (b) illustrates the normalized MAEs of water saturation at 5 cm beneath surface, as well 559 as temperature at surface and 1 m depth. The assumption of equal ice-liquid density primarily 560 affects the estimation of water saturation profile. It can lead to about 5% ~ 15% error relative to 561 the annual change range, and the error percentage tends to slightly decrease when applying 562 hillslope-scale models due to the inclusion of lateral flow. Apart from this, neglecting soil 563 cryosuction still makes the largest impact. Surface temperature is the least affected variable among 564 all these model outputs even if cryosuction is not included in modeling. However, at 1 m depth, 565 error can increase to 10% ~ 30% by simulation without cryosuction representation.

Conclusion 574
Simplification of Arctic process representation is an essential consideration when developing 575 process-rich models for thermal permafrost hydrology. There are three subsurface processes that 576 are commonly described in a simplified approach for many Arctic tundra models: (i) ice is 577 prescribed the same density as liquid water; (ii) the effect of soil cryosuction is neglected; (iii) 578 advective heat transport is neglected. Here we investigated the influence of these simplified 579 representations on modeling field-scale permafrost hydrology. 580 To do this, we conducted an ensemble of simulations using the Advanced Terrestrial Simulator 581 (ATS v1.2) to evaluate the impact of the above three process simplifications on field-scale 582 predictions. The ensemble of simulations consisted of 62 numerical experiments considering 583 various conditions, including different climate conditions and soil properties at three sites of 584 Alaska, and different model scale conceptualizations. For evaluation, we compared integrated 585 variables (evaporation, discharge), averaged thaw depth, and pointwise variables (temperature, 586 total water saturation), among different models to access the deviation of applying a simplified 587 modeling assumption. The main conclusions of this study are summarized as follows: 588 1) Excluding soil cryosuction in permafrost models can cause significant bias in most 589 hydrological variables. Especially, according to this study, the average deviation in 590 evaporation, discharge, and thaw depth may reach 10% ~ 20%, 50% ~ 60%, and 10% ~ 591 30%, respectively, relative to the corresponding annual average values. The prediction 592 error for discharge may grow if rainfall rates increase. In the case of pointwise variables, 593 the error in temperature increases from a small amount at the surface up to 10% ~ 30% at 594 1 m beneath surface. The prediction of subsurface temperature and water saturation is 595 especially affected when considering hillslope scale models. Therefore, soil cryosuction 596 should be included when modeling permafrost change. 597 2) Assuming equal ice-liquid density will not result in especially large deviations when 598 predicting most of the hydrological variables, particularly at hillslope scales. It primarily 599 affects the prediction of soil water saturation profile and can cause 5% ~ 15% error relative 600 to the annual saturation fluctuation range. This difference may have consequences for the 601 carbon cycle with regards to the production of methane versus carbon dioxide. Assigning 602 liquid water density for ice may reduce computational time to a small extent in ATS, 603 dependent on simulating conditions and spatial and temporal scales. 604 3) For a general Arctic tundra system, the prediction error in most variables after neglecting 605