Three-Dimensional Stefan Equation for Thermokarst Lake and Talik 1 Geometry Characterization 2

Abstract. Thermokarst lake dynamics, which plays an essential role in carbon release due to permafrost thaw, is affected by various geomorphological processes. In this study, we derive a three-dimensional (3D) Stefan equation to characterize talik geometry under a hypothetical thermokarst lake in the continuous permafrost region. Using the Euler equation in the calculus of variations, the lower bounds of the talik were determined as an extremum of the functional describing the phase boundary area with a fixed total talik volume. We demonstrate that the semi-ellipsoid geometry of the talik is optimal for minimizing the total permafrost thaw under the lake for a given annual heat supply. The model predicting ellipsoidal talik geometry was verified by talik thickness observations using transient electromagnetic (TEM) soundings in Peatball Lake on the Arctic Coastal Plain (ACP) of Alaska. The lake width-depth ratio of the elliptic talik can characterize the energy flux anisotropy in the permafrost although the lake bathymetry cross section may not be elliptic due to the presence of near-surface ice-rich permafrost. This theory suggests that talik development stabilizes thermokarst lakes by ground subsidence due to permafrost thaw while wind-induced waves and currents are likely responsible for the elongation and orientation of thermokarst lakes in certain regions such as the ACP of northern Alaska.


(13) 184 This expression indicates that the heat required for lake expansion is proportional to the weighted 185 phase boundary area. 186 2.2 Optimum phase boundary shape as extremum 187 The calculus of variation, often referred to as a functional analysis, is the mathematical technique 188 to find an extremum (minimum or maximum) of the system in terms of a function type instead of 189 a variable (e.g. Courant and Hilbert, 1954;Gelfand and Fomin, 1963). Here, we present the 190 thermally optimum function type ( , ) of the phase boundary using this method. As presented 191 in the previous section, the heat consumption rate for talik expansion is represented by the 192 weighted phase boundary area while the time-integrated heat supply is equivalent to the thawed 193 permafrost volume. Assuming heat thaws the most susceptible region of the permafrost near the 194 heat source first, the shape of a talik may minimize the total permafrost thaw with a given 195 amount of incoming energy. Hence, this variational principle states that the optimum talik shape 196 should minimize the phase boundary area for the total talik expansion. The weighted phase 197 boundary area A and its volume V can be expressed as follows: If the permafrost degradation rate is uniform and constant throughout the basin (Panel A: 239 uniform permafrost in Figure 2), the lake bathymetry tends to be an ellipsoid shape. However, as 240 the ice-rich layer (ice wedges) is typically developed near the surface on the ACP (e.g. 241 Kanevskiy et al., 2013Kanevskiy et al., , 2016, the bathymetry may have a flatter bottom like a rectangular cross 242 section (Panel B: layered permafrost in Figure 2) because the ice-rich layer is characterized by 243 much higher thaw settlement than the ice-poor permafrost at depth. Therefore, proportionality 244 between talik thickness and lake water depth or uniform permafrost is unlikely reasonable 245 assumption due to the ice rich layer presence. Indeed, Hinkel et al. (2012) showed many flat-246 bottomed lakes through the extensive bathymetry surveys across the ACP of Alaska using a 247 GPS-enabled sonar from a boat. 248 Additionally, as hydrology also affects the lake water level, the apparent lake bathymetry or lake 249 water depth, ℎ( , ) must be adjusted by the water loss (or gain) per unit area. Therefore, 250 where (m) is the elevation difference between the current water surface and original ground 252 surface before lake formation. At the lake center, 253 for the demonstrative model application in this study as it has been relatively well-documented. 260 Figure 3 shows the location of Peatball Lake as well as the subregions that will be presented later.  Figure 4 shows the observed talik thicknesses by the TEM 285 sounding (dots) and the fitted theoretical talik thickness estimates (contour lines) superimposed 286 over the corresponding lake bathymetry measured by Lenz et al., 2016. 287 The geometric parameters of the semi-ellipsoid model such as the talik center depth (D), the 288 cross-sectional aspect ratios (αx & αy), lake orientation azimuth angle and the lake center 289 location were systematically determined by grid searching to minimize the root mean square 290 difference (RMSD) between the model and thaw front obtained from the TEM data. The 291 optimum parameters for the smallest RMSD (5.94 m) are shown in Figure 4. Unexpectedly, the 292 basin orientation angle was found to be 23 degrees east from true north, unlike the orientation of 293 other surrounding lakes in the region. Comparison between the extrapolated talik geometry and 294 the lake bathymetry (Lenz et al., 2016) suggests the possibility of coalescence of two basins in 295 the past; a relatively common occurrence on the ACP of Alaska. However, if we had more TEM 296 measurement points, particularly in the "possible talik sub-basin", the fitted talik geometry could 297 be different as the model was only fitted for the 27 TEM soundings. Lake taliks tend to have a 298 semi-ellipsoidal shape, at least locally, as indicated by the very good fit of the elliptic model to 299 the TEM measured talik thicknesses (see Figure 4 with overall RMSD = 5.94 m). The idealized, thermally optimum geometry model can be used to analyze lake formation history of the 301 irregular talik associated with multi-generation lakes such as Peatball Lake. 302 Additionally, the gaps between the shoreline and the modeled talik extent located along north 303 and east shores occur where lake expansion is occurring most rapidly (Lenz et al., 2016). It has 304 been reported that thaw lake banks continuously retreat through a combination of thermal and 305 mechanical processes, although there is significant variability in rate of bank retreat depending 306 on region (Hopkins, 1949;Hopkins et al., 1955;Tomirdiaro, 1982;Rampton, 1982 and seasonal snow cover (Ling and Zhang, 2003a). The disagreement between the lake and talik 310 extents on the north and east shores of Peatball lake implies that rapid horizontal lake expansion 311 can locally dominate permafrost thaw and subsidence processes even in a lake with a talik. 312 Figure 5 compares the observed lakebed and talik profiles in Peatball Lake along the north-south 313 center line and along transects (b) -(c), respectively. Note that the TEM transects for the talik is 314 not a straight line (See figure 4); therefore, the fitted theoretical line shows irregularity. Figure  315 5A illustrates that the lakebed profile is characterized by flatter trapezoidal geometry compared 316 to the elliptic talik. In fact, there is a clear inflection in the linear regression line at a talik depth 317 of ~50 m in Figure 5B. From the slopes of the regression lines, the permafrost degradation rates 318 are computed as 97.3 % and 99.7 % for the shallow talik section (50 m or less) and the deep 319 section (50 m or more), respectively. This analysis suggests that the subsidence due to 320 permafrost thaw continues even after the shallow ice-rich part of the permafrost (about 4 meters, 321 Kanevskiy, 2011) is thawed while it has diminished around the depth of 50 meter under Peatball 322 lake. This case study demonstrates a link between lake bathymetry and talik thickness associated 323 with a layered permafrost structure and the wind erosion effect. This simple expression may be a useful tool to link the lake depth-width ratio, the lake average 345 heat flux , and the vertical temperature gradient at the base of the talik. Since < 0 in 346 the permafrost near the talik boundary, the D/a is always less than 1 (flatter than a semi-sphere). 347 However, the depth-width ratio of the talik depends on the vertical temperature slope near the 348 talik boundary, which is likely affected by talik age. For instance, Mackay's analytical model 349 (1962) suggests that the vertical temperature gradient below the lake center begins steeply at the 350 talik initiation, and then over time it approaches a lower slope at equilibrium. Therefore, the 351 formula in Equation (29) suggests that a younger talik should be flatter while an older talik 352 approaches a deeper semi-spheric shape (D/a → 1). 353 Table 1 shows the estimated incoming heat flux with the key parameters using the proposed 354 formula (Equation (29)

373
To illustrate the applicability of the thermal model presented here, the available hypothetical 374 models of thermokarst lake growth are compiled in Figure 6. This diagram focuses on the 375 physical processes after the lake initiation stage assuming the bio-ecological effects are 376 negligible. 377 Figure 6 illustrates the evolution of the talik in ice-rich permafrost over time, with driving 378 processes shown in the right panel. In Stage A, the mechanical processes of wave erosion and 379 thaw slumping along lake margins dominate lake expansion in summer, and shallow water favors 380 grounded lake ice in the winter. In time (Stage B), the lake deepens from thaw subsidence 381 beneath the older lake center. Winter ice may freeze to the lake bed, but heat loss is insufficient 382 to freeze the underlying thawed lake bed sediments. A shallow talik develops as thermal 383 processes work in tandem with mechanical processes, the latter now enhanced by more vigorous 384 lake circulation. By Stage C, the talik is well developed beneath the entire lake basin as ground 385 subsidence continues. Eventually (Stage D), the winter ice cover no longer extends to the lake 386 bed, but instead floats atop a residual pool of lake water, while milder vertical temperature 387 gradient beneath the lake deepens the talik as the lake matures. Thermomechanical erosion of 388 lake margins, especially if there are prominent banks in hilly terrain, promotes sedimentation on 389 near-shore shelves, and the underlying talik may begin refreezing. If the lake hasn't drained by 390 this point (Stage E), the talik beneath the lake center extends deeper into the permafrost although 391 subsidence may cease as the excess ice content diminishes with depth. Where many large, old 392 lakes exist, the permafrost may be riddled with deep taliks, and some may eventually penetrate to 393 the permafrost base to create a through-going talik. 394 Talik development is a natural processes governed by local conditions, but can be impacted by 395 direct and indirect human activities. Conditions that favor talik initiation and growth include: 396 • Deepening lake waters triggered by greater precipitation and/or reduced evaporation, 397 which promotes a floating ice regime 398 • Presence of ice-rich sediments (e.g., Yedoma) beneath lakes 399 • Warmer lake water induced by regional warming or by longer ice-free summers 400 • Thinner winter ice cover due to warmer winter temperatures and/or deeper snow 401 Conversely, talik growth cessation or contraction can occur when these same drivers are reversed, 402 if the lake partially or completely drains, or when the lake basin is filled with sediments. The 403 latter scenario is more likely in hilly terrain when the expanding lake erodes high banks and lake 404 currents redistribute sediments. The analytical expression of the lake geometry may be useful to analyze the horizontally oriented 408 lakes as well. From Equation (27) and (1) heat conduction from the lake water body. Thus, the lake aspect ratio may be written as, 418

Incoming radiation imbalance effect 420
One of the incoming surface energy flux inequalities may be caused by shortwave radiation 421 along the lake shoreline. The daily potential solar irradiation on a sloping surface can be 422 computed by the trigonometric function (e.g. Equation B.11 in DeWalle and Rango, 2008). The 423 total daily radiation is a function of latitude and bank slope angle, which depends on the Figure 7 shows the computed mean daily potential solar irradiation on the sloping lakeshore (I'q) 426 relative to a flat surface (Iq) during the summer period (June-August) at three different latitudes. 427 The shape of this diagram may correspond to the shape of a thermokarst lake as the enhanced 428 radiation results in more permafrost thaw. The difference in relative incoming radiation will 429 diminish as bank slope angle lessens. In general, the south facing slope along the northern shore 430 tends to receive more radiation than the north facing slope (e.g. Séjourné et al., 2015). This 431 tendency is more pronounced in lower latitude zones due to the higher mid-day sun angle. 432 It is interesting that at 65 and 60 degree latitude the north and south facing banks receive slightly 433 less radiation than east and west facing slopes, while an opposite result occurs at 70 degree 434 latitude (Figure 7). Therefore, the radiation imbalance may partially explain the north-south  Figure 19). However, 437 because these small differences in incoming radiation imbalance alone are insufficient to result 438 in the distinctive lake elongation in the ACP, they likely introduce rather minor additional 439 complexities in lake spatial shape. The shallow wave theory states that the sediment mobilization due to wind wave only occurs in 475 shallow water (wave height >4 % of water depth, e.g. Reeve et al., 2018). Therefore, the 476 contribution of the wind wave effect to lake elongation may be reduced as the lake deepens. 477 Figure 9 shows a plot of lake length:width ratios versus the percent of lakes with a bedfast ice 478 regime in seven study regions in Alaska determined with satellite-based synthetic aperture radar 479 imagery (Engram et al., 2018). The study regions represent differences in permafrost 480 characteristics and climate that appear to be reflected in this comparison of length:width ratio 481 and the percent of lakes in a region that freeze to their bed and thus likely do not have a sub-lake 482 talik. For example, lakes in the Teshekpuk Lake and Kuparuk study areas have a shape that is 483 nearly twice as long as it is wide. In both of these regions, more than 80% of the lakes freeze to 484 their bed and likely do not have a talik. This is contrast to lakes located near Umiat and on the 485 Seward Peninsula, that have primarily developed in Yedoma permafrost deposits. Lakes near 486 Umiat and on the Seward Peninsula tend to be more circular (L:W = 1.3 to 1.4) and more than 487 90% likely have a talik since they do not freeze to their bed in the winter. The differences 488 observed here relative to elongation of lakes and whether the region primarily has lakes that 489 freeze to their bed or not likely demonstrates a key aspect related to the role of wind-wave 490 erosion. In general, the shallower lakes common in coastal areas, such as Teshekpuk, Barrow,491 and Kuparuk, are more elongated likely due to wind wave erosion. Whereas lakes in Umiat and 492 Inigok with a thicker ice-rich permafrost layer tends to rounder because of more rapid talik 493 development and more subsequent subsidence. This remote-sensing based evidence implies that 494 the wind effect seems to be limited by the lake thermal subsidence due to talik development 495 underneath while the lake with the bedfast ice may continue elongating by the wind erosion. be the result of topography (e.g. slope and aspect of the ground surface) while lake elongation is 503 likely caused by wind wave erosion. As described above, preferential bank thaw at the lake ends 504 can be explained by the insulation effect of the sediments carried by the water current (likely, the 505 C&H type circulation) because the sublittoral shelf may be initiated at this stage. 506 When seasonal thawing penetrates more deeply than the annual freezing depth, a talik may be 507 initiated, typically around the deepest point near the center of the lake (Lachenbruch et al. 1962). 508 Sellmann (1975) described this process, which is one of the mechanisms for shelf formation in a 509 thermokarst lake. For the horizontal expansion stage A in Figure 6, the proposed quasi-steady 510 state thermal model may not be appropriate because the lakeshore expansion imbalance occurs at 511 least minimally throughout the lake expansion process. However, the 3D Stefan equation may 512 be able to characterize the talik in the initiation stage B in Figure 6. 513 Once the talik is established, the 3D Stefan's thermal model proposed here suggests that the talik 514 may begin to influence lake geometry. Since sediment mobilization due to wind-driven waves 515 occurs in shallow water, lake elongation by waves may diminish as the lake deepens via ground 516 subsidence (Figure 9). Lake water effectively collects energy from the surface during summer 517 and the talik stores the excess heat throughout the winter. occurs before talik formation, the horizontal lake characterization derived in this study may not 522 be fully applicable to the analysis of thaw lakes on the ACP. In fact, the disagreement of the 523 talik and lake extents in Peatball Lake application illustrates the multiple effects on the lake 524 bathymetry and orientation. Clearly, however, talik expansion and concurrent subsidence 525 stabilizes lake geometry and contributes to lake roundness. 526 The applicability of the proposed 3D Stefan equation must be limited for lakes with high 527 sediment influx and for lakes with through talik. The paired sublittoral shelves on both lake sides 528 are commonly found in the sand dune areas of the southern ACP. The talik shape is likely altered 529 by uneven sediment deposition that affects the temperature gradient normal to the phase 530 boundary as mentioned by several researchers (Mackay, 1992; Hunter et al., 1981; West and 531 Plug, 2008). The shelves created by sediment redistribution due to lake water circulation adds 532 complexity to the ellipsoidal talik shape described in this study. Finally, if the talik penetrates 533 through the permafrost and becomes a throughgoing talik (Hinkel and Arp, 2015), the proposed 534 thermal theory herein is no longer applicable for thermokarst lake and talik characterization. 535

536
The theory presented here addresses the origin of the thermokarst lake ellipticity on the ACP. 537 Elliptic lake geometry results from minimizing overall thawing energy consumption for a given 538 incoming energy load. This is particularly applicable for mature, deep thermokarst lakes with 539 well-developed taliks. Additionally, existing hypothetical models were reviewed to illuminate 540 the thermal effect on the thermokarst lake morphology. 541 The derived ellipsoid talik model integrates the atmospheric forcing (or incoming energy), the 542 vertical thermal gradient, the thermal diffusivity of the permafrost, and the talik geometry. Heat 543 flux by conduction into the permafrost depends on the heat gradient of the underlying permafrost 544 (Fourier's law). As the vertical temperature slope diminishes with talik maturation, the depth-545 width ratio of the talik becomes larger creating a deeper talik; thus, much of incoming energy is 546 likely consumed for vertical rather than horizontal expansion. Conversely, during the early stages, 547 thermo-mechanical processes such as wind-driven wave erosion dominates horizontal expansion 548 and elongation of the lake. Consequently, this theory elucidates how talik expansion and 549 concurrent permafrost degradation stabilizes the shape of thermokarst lake to one that is more 550 round rather than elliptical. 551 The semi-ellipsoidal 3D Stefan equation is, to our knowledge, the first geometric model 552 explicitly derived only from the energy conservation equation at the phase boundary. The vector 553 form of the energy conservation equation (Equation 5) in a 3D anisotropic thermal field was 554 integrated at the phase boundary area under the isolated general-shaped lake to quantify the total 555 energy balance. It was shown that the total lake fusion energy or lake expansion rate is 556 equivalent to the weighted phase boundary area. The optimum talik shape function was 557 determined by the variational principle as an extremum of the functional that minimizes the total 558 thawing energy consumption. Thus, the resultant semi-ellipsoid equation (Equation 22) can be 559 considered the 3D Stefan equation because it describes the optimum geometry of phase boundary. 560 The derived semi-ellipsoid function was applied to Peatball Lake, ACP of Alaska, where the 561 talik was extensively surveyed using TEM soundings. The pure geometric fitting exercise met 562 the 27 measured TEM data point well with RMSD of 5.94 m, although the talik orientation 563 disagreed with orientation of Peatball Lake and other surrounding lakes. This may be induced by 564 the irregularity due to the rapid and uneven horizontal lake expansion, or possibly by basin 565 coalescence. Comparing the observed talik thickness to the observed lake bathymetry indicated 566 two distinctive permafrost degradation scenarios: significant subsidence by near-surface ice-rich 567 layer thaw and minor contribution of subsidence due to ice-poor permafrost thaw at depth. 568 Consequently, lake water depth is affected by uneven subsidence of thawing permafrost, the 569 interannual water balance; the spatial lake shape irregularity was determined during earlier stage 570 of development. Therefore, careful consideration is required for the analysis of the relationship 571 between lake bathymetry and talik thickness. Nevertheless, this theoretical technique can be 572 used as guidance to partition various effects such as talik development and thaw subsidence, 573 wind wave erosion, lake ice thickness, surficial geology type, and sediment transport by lake 574 water current. 575 Appendix A: Alternative derivation using isoperimetric inequality Eliminating D from these expressions yields, 624 where M is a positive constant. Therefore, as 626

887
Lakes that are more circular in shape tend to occur where the majority of lakes in an area do not freeze to 888 their bed and thus likely have a sub-lake talik.