An improved laminar flow inversion method for ice thickness based on glacier bottom sliding

Abstract

This invention discloses an improved laminar flow inversion method for ice thickness based on glacier bottom sliding. The method prepares glacier surface velocity and digital elevation model (DEM) data in raster format, and glacier boundary data in vector format. The glacier surface slope is calculated based on the DEM, and a glacier mask is generated based on the glacier boundary data. These five data types are then unified into a projected coordinate system. A laminar flow theory formula is constructed for calculating glacier thickness H. The initial glacier thickness is obtained by assuming the glacier bottom sliding velocity is zero. The valley shape factor, bottom shear stress, and bottom sliding velocity are sequentially updated to obtain a new glacier thickness. The method calculates the difference between adjacent glacier thicknesses. When the difference meets a set condition, the calculation is terminated and the result is output. Otherwise, the difference is assigned to the initial glacier thickness and the calculation is restarted until the difference meets the set condition. This method effectively improves the accuracy of ice thickness models based on traditional laminar flow theory and can also automatically obtain subglacial topography.

Description

Technical Field

[0001] The present invention belongs to the research field of mountain glacier thickness estimation methods, and in particular relates to an improved laminar flow inversion method for ice thickness based on glacier bottom sliding. Background Art

[0002] Glaciers are natural, long-standing ice masses in polar and alpine regions, formed by the repeated processes of compaction, crystallization, and freezing of solid precipitation, such as snow and hail. Climatic conditions such as precipitation and temperature, as well as topographical conditions at the altitude of mountains above the snowline, jointly determine the formation of mountain glaciers, directly influencing their number and size. Globally, glaciers cover approximately 16 million square kilometers, accounting for approximately 11% of the total land area and consolidating 68.7% of the world's freshwater resources. Mountain glaciers are primarily found in the Alps, Alaska, high mountain Asia, and the southern Andes, flowing slowly downward along valleys from high to low. As a crucial component of the cryosphere, glaciers are not only an indicator of climate change but also a key constraint on human socioeconomic development. Due to continued global warming in recent decades, glaciers in alpine Asia have experienced significant decline and mass loss, contributing to accelerated sea level rise. Furthermore, rapid glacier movement can trigger various glacial hazards, such as glacial lake outburst floods (GLOFs) and ice avalanches. Therefore, developing ice thickness inversion models to clarify the thickness distribution and volume of mountain glaciers is a crucial foundation for glaciological research.

[0003] Currently, regional and global-scale mountain glacier thickness and volume estimates primarily utilize glacier topography, kinematic characteristics, and geometric features, inverted using glacier dynamics. However, significant variations exist in the inversion of ice volume using existing laminar flow theory-based ice thickness models. These variations are attributed to the accuracy of the ice thickness models, the quality of the input data, and parameter estimation. In these models, laminar flow models estimate ice thickness based on the glacier's surface velocity and basal slip. Typically, laminar flow models assume a scaling relationship between basal slip and surface velocity when estimating ice thickness. For example, a basal slip-to-surface velocity ratio of 25% (i.e., the ratio of slope to surface velocity) or zero basal slip in winter can be assumed. However, basal slip is affected by factors such as subglacial topography, basal shear stress, and ice cover. The ratio of basal slip to surface velocity varies among glaciers, with measured values ​​ranging from 0.03 to 1. Therefore, incorporating basal slip into laminar flow theory is crucial for improving the accuracy of ice thickness estimates. Glacier surface velocity can be considered the sum of basal slip and ice deformation. Glacier basal slip is closely related to the base temperature and water content at the ice-rock interface. Subglacial terrain is generally irregular and has numerous obstacles. Glacier movement is driven by two mechanisms: pressure melting and creep enhancement. However, obtaining information about ice bed geometry and temperature remains challenging. Due to the technical complexity of field measurements, slip laws are commonly used to model basal sliding in glacier dynamics studies. In 1957, Veltman developed a slip law that assumes the glacier rests on solid bedrock, with no shear stress or cavities at the interface, and the ice bed is idealized as a cubic convexity. However, basal shear stress, effective pressure (the difference between the water pressure of the overlying layer and the base of the ice), and ice bed roughness all affect basal sliding. Therefore, the slip law needs to account for bed cavities and incorporate constraints to prevent the ratio of basal shear stress to effective stress from increasing with increasing sliding velocity. Coulomb's friction law satisfies these conditions, causing the ratio of basal shear stress to effective stress to initially increase to a maximum before gradually decreasing and stabilizing. Therefore, the slip law has been used to determine basal sliding in laminar flow models to improve the accuracy of ice thickness estimates, but further research is lacking. Summary of the Invention

[0004] This invention addresses the problem of low ice thickness estimation accuracy caused by the failure of existing laminar flow-based ice thickness models to account for bottom sliding. By proposing an improved laminar flow-based ice thickness inversion method based on glacier bottom sliding, this method is implemented by incorporating the law of glacier bottom sliding. This method enables spatial estimation of glacier bottom sliding, and further calculation of glacier thickness using laminar flow theory. This helps improve ice thickness inversion accuracy and also enables automated acquisition of subglacial topography, providing an important data foundation for glacier dynamics research.

[0005] The above-mentioned purpose of the present invention is achieved by the following technical means:

[0006] An improved laminar flow inversion method for ice thickness based on glacier bottom sliding is characterized by comprising the following steps:

[0007] Step 1: Obtain the glacier surface velocity u s , digital elevation model DEM, glacier boundary data, glacier surface slope α, glacier mask and unified into a projected coordinate system;

[0008] Step 2: Construct a model containing the glacier surface velocity u s and the sliding speed u at the bottom of the glacier b The laminar flow theory formula is used to construct the calculation formula for the glacier thickness H;

[0009] Step 3: Sliding speed u at the bottom of the glacier b Assign a value of 0m / s, and use the glacier thickness calculation formula constructed in step 2 to obtain the initial glacier thickness H0. At the same time, use the digital elevation model (DEM) to subtract the initial glacier thickness H0 to obtain the initial subglacial terrain Z0.

[0010] Step 4: Update the valley shape factor f based on the average glacier width w and glacier thickness H;

[0011] Step 5: Update the bottom shear stress τ using the updated valley shape factor f b ;

[0012] Step 6: Based on the updated bottom shear stress τ b Sliding speed u at the bottom of the glacier b Make updates;

[0013] Step 7: Based on the valley shape factor f updated in step 4 and the glacier bottom sliding velocity u updated in step 6 b , using the calculation formula of glacier thickness to get the new glacier thickness H i ;

[0014] Step 8: Calculate the difference in glacier thickness between two adjacent calculations ΔH = H i -H i-1 , i≥1, i is the number of times the glacier thickness is calculated. When the difference ΔH is less than the set threshold, the last calculated glacier thickness is output, and the final subglacial terrain is obtained by subtracting the last calculated glacier thickness from the digital elevation model DEM. The last updated glacier bottom sliding speed is used as the final glacier bottom sliding speed; when the difference ΔH is greater than or equal to the set threshold, the glacier thickness H is used. i Take H0 as the new initial glacier thickness and return to step 3.

[0015] In step 1, the glacier surface slope α is obtained by calculating the digital elevation model DEM, and the glacier mask is obtained by the glacier boundary data. The digital elevation model DEM, glacier boundary data, glacier surface slope α, and glacier mask are resampled to the glacier surface velocity u s Same resolution.

[0016] In step 2, the laminar flow theory formula is based on the following formula:

[0017]

[0018] Where: H represents the thickness of the glacier, u s and u b are the surface velocity and bottom sliding velocity of the glacier, n is the Glen flow law index; A is the ice creep coefficient; τ b represents the shear stress at the base of the glacier;

[0019] Glacier bottom shear stress τ b Based on the following formula:

[0020] τ b =fρgHsinα

[0021] Where f is the valley shape factor, ρ is the density of the ice body, g is the acceleration of gravity, and α is the slope of the glacier surface;

[0022] The shear stress at the bottom of the glacier τ b Substituting the calculation formula into the laminar flow theory formula, we get the calculation formula for glacier thickness:

[0023]

[0024] The valley shape factor f is updated in step 4 based on the following formula:

[0025]

[0026] Where w is the average width of the glacier.

[0027] The average width w of the glacier is obtained based on the following steps:

[0028] The points on the glacier centerline are extracted from the glacier mask, and the coordinates of each point are calculated. Then, the twice the distance from each point to the glacier boundary in the direction perpendicular to the glacier centerline is calculated as the glacier width, and the average of the glacier thicknesses obtained by each calculation is taken to obtain the average width w of the glacier.

[0029] The sliding speed u of the glacier bottom in step 6 is b Update based on the following formula:

[0030]

[0031] Where N is the effective pressure of the glacier, C is the maximum ratio of the shear stress at the bottom of the glacier to the effective pressure, and A is the maximum value of the ratio of the shear stress at the bottom of the glacier to the effective pressure. s is the glacier sliding parameter, B is the mobility parameter, n is the Glen flow law index, λ is the wavelength of the glacier topography; r is the roughness of the glacier.

[0032] In step 6, the effective pressure N of the glacier is based on the following formula:

[0033] N=ρgH·cosα

[0034] Where ρ represents the density of the ice body, g is the acceleration of gravity, H is the thickness of the glacier, and α is the slope of the glacier surface;

[0035] The maximum value of the ratio of shear stress to effective pressure at the bottom of the glacier, C, is based on the following formula:

[0036] C=0.84±0.02m max

[0037] Among them, m max To obtain the maximum slope of the subglacial terrain based on the initial subglacial terrain Z0;

[0038] The glacier sliding parameter is based on the following formula:

[0039]

[0040] Wherein, B is the fluidity parameter, B = 430 ± 40 MPa -3 yr -1 , λ is the wavelength of the ice bed topography, r is the roughness of the ice bed;

[0041] The roughness of a glacier is based on the following formula:

[0042]

[0043] Among them, k represents the number of pixels of the ice bed, j represents the sequence number of the pixel, and z j represents the elevation of the jth pixel of the ice bed, Indicates the average elevation of the ice bed.

[0044] A computer device includes a memory and a processor, wherein the memory stores a computer program, and the processor implements each step of the above-mentioned ice thickness inversion method when executing the computer program.

[0045] A computer-readable storage medium stores a computer program, which implements the various steps of the above-mentioned ice thickness inversion method when executed by a processor.

[0046] A computer program product includes a computer program, which implements the various steps of the above-mentioned ice thickness inversion method when executed by a processor.

[0047] Compared with the prior art, the present invention has the following beneficial effects:

[0048] (1) High precision: The present invention uses glacier surface velocity, surface elevation, and glacier boundary data as input, determines bottom sliding through the initial glacier thickness and initial subglacial topography estimated by traditional laminar flow theory, and inverts ice thickness through cyclic iteration to laminar flow theory, effectively improving the accuracy of the ice thickness model based on traditional laminar flow theory.

[0049] (2) Practicality: The present invention can invert the glacier thickness to obtain high-precision data and can also obtain the subglacial topography, providing a basis for the dynamic study of mountain glaciers. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] Figure 1 It is a schematic diagram of the process of the present invention;

[0051] Figure 2 is the glacier surface velocity of ChhotaShigri glacier;

[0052] Figure 3 The glacier base sliding velocity of ChhotaShigri glacier estimated by this method;

[0053] Figure 4 is the ratio of the glacier base sliding velocity to the surface velocity of ChhotaShigri glacier;

[0054] Figure 5 The thickness distribution of ChhotaShigri glacier estimated by this method;

[0055] Figure 6 Subglacial topography of the ChhotaShigri glacier estimated using this method. DETAILED DESCRIPTION

[0056] In order to facilitate those skilled in the art to understand and implement the present invention, the present invention is further described in detail below in conjunction with embodiments. It should be understood that the embodiments described herein are only used to illustrate and explain the present invention and are not used to limit the present invention.

[0057] Example 1:

[0058] like Figure 1 As shown, a method for inverting ice thickness based on improved laminar flow at the bottom of a glacier includes the following steps:

[0059] Step 1: Obtain the glacier surface velocity u in raster format sThe glacier boundary data in vector format is obtained by using the digital elevation model (DEM). The glacier surface slope α is calculated based on the digital elevation model DEM. The glacier mask is generated based on the glacier boundary data. The glacier surface velocity u is converted to s , digital elevation model, glacier boundary data, glacier surface slope, and glacier mask are unified into a projected coordinate system; to ensure the consistency of pixel size, the digital elevation model DEM, glacier boundary data, glacier surface slope α, and glacier mask are resampled to the glacier surface velocity u s Same resolution.

[0060] In this embodiment, the ChhotaShigri glacier is used as the object for calculation. Figure 2 is the glacier surface velocity u of ChhotaShigri Glacier s .

[0061] Step 2: Based on the glacier surface velocity u s and the sliding speed u at the bottom of the glacier b , construct the laminar flow theory formula:

[0062]

[0063] Where: H (m) represents the thickness of the glacier; u s and u b are the surface velocity and bottom sliding velocity of the glacier, respectively; n is the Glen flow law index, which is 3; A is the ice creep coefficient, which is related to the ice temperature, structure, and water content, and is 2.4×10 -24 Pa -3 s -1 ; τ b represents the shear stress at the base of the glacier, calculated using the shallow ice approximation formula.

[0064] τ b =fρgHsinα (2)

[0065] Where: f represents the valley shape factor, which is related to the glacier cross section, driving stress, and bottom stress, and its initial value is 0.8; ρ represents the density of the ice body, which is 900 kg·m -3 ; g is the acceleration due to gravity, which is 9.8 m·s -2 α is the slope of the glacier surface. Substituting equation (2) into equation (1) yields the final formula for calculating the glacier thickness H.

[0066]

[0067] Step 3: First, slide at the bottom of the glacier at a speed u b If the sliding speed at the bottom of the glacier is unknown, then bThe initial glacier thickness H0 is obtained by using the calculation formula for glacier thickness H constructed in step 2. At the same time, the initial subglacial terrain Z0 is obtained by subtracting the initial glacier thickness H0 from the digital elevation model (DEM).

[0068] Step 4: Update the valley shape factor f based on the average glacier width w and glacier thickness H. The valley shape factor f is updated and calculated based on the glacier mask, glacier thickness, and glacier boundary data.

[0069]

[0070] Where w is the average width of the glacier.

[0071] The average width w of the glacier is obtained as follows: points on the glacier centerline are extracted based on the glacier mask, and the coordinates of each point are calculated; then, the glacier width is calculated as twice the distance from each point to the glacier boundary in the direction perpendicular to the glacier centerline, and the average value is taken to obtain the average width w of the glacier.

[0072] Step 5: Update the bottom shear stress τ according to the updated valley shape factor f b , calculated according to the shallow ice approximate formula of formula (2).

[0073] Step 6: According to the bottom sliding law of the glacier, based on the updated bottom shear stress τ b Sliding speed u at the bottom of the glacier b to update.

[0074]

[0075] Where N represents the effective pressure of the glacier, which is the difference between the ice load and the water pressure at the bottom of the glacier. When the bottom water pressure is ignored, the calculation formula is N = ρgH·cosα; C is the maximum value that can be achieved by the ratio of shear stress to effective pressure at the bottom of the glacier, C = 0.84 ± 0.02m max , where m max Represents the maximum slope of the subglacial terrain. The maximum slope m of the subglacial terrain is obtained based on the initial subglacial terrain Z0. max ; A s is the glacier sliding parameter, which is affected by the roughness of the subglacial terrain and the temperature of the ice body. Where B is the fluidity parameter. Assuming that the ice is isothermal, B = 430 ± 40 MPa -3 yr -1 ; n is the Glen flow law index, which is 3; λ is the wavelength of the ice bed topography, r is the roughness of the ice bed, which is calculated using the root mean square height (RMSH), Where k is the number of pixels on the ice bed, j is the pixel number, and z is the pixel number. j represents the elevation of the jth pixel of the ice bed, Represents the average elevation of the ice bed; after calculating various parameters, we put them into formula (5) to obtain the updated sliding speed u at the bottom of the glacier b .

[0076] Figure 3 The sliding velocity u of the bottom of the ChhotaShigri glacier estimated by this method is b .

[0077] Step 7: Update the sliding speed u of the glacier bottom updated in step 6 b , the valley shape factor f updated in step 4 is introduced into the glacier thickness calculation formula in step 2 to obtain the new glacier thickness H i .

[0078] Step 8: Calculate the difference in glacier thickness between two adjacent calculations ΔH = H i -H i-1 , i≥1, i is the number of times the glacier thickness is calculated. When ΔH is less than the set threshold, it means that the condition is met, the calculation is terminated, and the last calculated glacier thickness is output. The final subglacial terrain is obtained by subtracting the last calculated glacier thickness from the digital elevation model DEM, and the last updated glacier bottom sliding speed is used as the final glacier bottom sliding speed. When the difference ΔH is greater than or equal to the set threshold, the glacier thickness H is used as the final glacier bottom sliding speed. i Let the value assigned to the new initial glacier thickness H0 and the glacier bottom sliding speed u updated according to step (6) b , return to step 3.

[0079] Figure 5 The thickness distribution of ChhotaShigri glacier calculated by this method, Figure 6 The final subglacial topography of the ChhotaShigri Glacier calculated using this method.

[0080] Those skilled in the art will understand that all or part of the processes in the above-mentioned embodiment methods can be implemented by instructing related hardware through a computer program. The computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above-mentioned methods.

[0081] Example 2:

[0082] In this embodiment, a computer device is further provided, including a memory and a processor. The memory stores a computer program, and the processor implements the steps in the above embodiment 1 when executing the computer program.

[0083] Example 3:

[0084] In this embodiment, a computer-readable storage medium is provided, on which a computer program is stored. When the computer program is executed by a processor, the steps in the above-mentioned embodiment 1 are implemented.

[0085] Example 4:

[0086] In this embodiment, a computer program product is provided, including a computer program. When the computer program is executed by a processor, the steps in the above-mentioned embodiment 1 are implemented.

[0087] It should be noted that the embodiments described herein are merely illustrative of the spirit of the present invention. Persons skilled in the art may make various modifications, additions, or substitutions to the described embodiments without departing from the spirit of the present invention or exceeding the scope of the appended claims.

Claims

1. A method for inverting ice thickness based on improved laminar flow at the bottom of glacier, characterized by: The following steps are involved: Step 1: Obtain the glacier surface velocity u s , digital elevation model DEM, glacier boundary data, glacier surface slope α, glacier mask and unified into a projected coordinate system; Step 2: Construct a model containing the glacier surface velocity u s and the sliding speed u at the bottom of the glacier b The laminar flow theory formula is used to construct the calculation formula for the glacier thickness H; Step 3: Sliding speed u at the bottom of the glacier b Assign a value of 0m / s, and use the glacier thickness calculation formula constructed in step 2 to obtain the initial glacier thickness H0. At the same time, use the digital elevation model (DEM) to subtract the initial glacier thickness H0 to obtain the initial subglacial terrain Z0. Step 4: Update the valley shape factor f based on the average glacier width w and glacier thickness H; Step 5: Update the bottom shear stress τ using the updated valley shape factor f b ; Step 6: Based on the updated bottom shear stress τ b Sliding speed u at the bottom of the glacier b Make updates; Step 7: Based on the valley shape factor f updated in step 4 and the glacier bottom sliding velocity u updated in step 6 b , using the calculation formula of glacier thickness to get the new glacier thickness H i ; Step 8: Calculate the difference in glacier thickness between two adjacent calculations ΔH = H i -H i-1 , i≥1, i is the number of times the glacier thickness is calculated. When the difference ΔH is less than the set threshold, the last calculated glacier thickness is output, and the final subglacial terrain is obtained by subtracting the last calculated glacier thickness from the digital elevation model DEM. The last updated glacier bottom sliding speed is used as the final glacier bottom sliding speed; when the difference ΔH is greater than or equal to the set threshold, the glacier thickness H is used. i As the new initial glacier thickness H0, return to step 3; In step 2, the laminar flow theory formula is based on the following formula: Where: H represents the thickness of the glacier, u s and u b are the surface velocity and bottom sliding velocity of the glacier, n is the Glen flow law index; A is the ice creep coefficient; τ b represents the shear stress at the base of the glacier; Glacier bottom shear stress τ b Based on the following formula: t b =fρgHsina Where f is the valley shape factor, ρ is the density of the ice body, g is the acceleration of gravity, and α is the slope of the glacier surface; The shear stress at the bottom of the glacier τ b Substituting the calculation formula into the laminar flow theory formula, we get the calculation formula for glacier thickness:

2. The ice thickness inversion method based on improved laminar flow of glacier bottom sliding according to claim 1 is characterized in that: In step 1, the glacier surface slope α is obtained by calculating the digital elevation model DEM, and the glacier mask is obtained by the glacier boundary data. The digital elevation model DEM, glacier boundary data, glacier surface slope α, and glacier mask are resampled to the glacier surface velocity u s Same resolution.

3. The ice thickness inversion method based on improved laminar flow of glacier bottom sliding according to claim 1 is characterized in that: The valley shape factor f is updated in step 4 based on the following formula: Where w is the average width of the glacier.

4. The ice thickness inversion method based on improved laminar flow of glacier bottom sliding according to claim 3 is characterized in that: The average width w of the glacier is obtained based on the following steps: The points on the glacier centerline are extracted from the glacier mask, and the coordinates of each point are calculated. Then, the twice the distance from each point to the glacier boundary in the direction perpendicular to the glacier centerline is calculated as the glacier width, and the average of the glacier thicknesses obtained by each calculation is taken to obtain the average width w of the glacier.

5. The ice thickness inversion method based on improved laminar flow of glacier bottom sliding according to claim 4 is characterized in that: The sliding speed u of the glacier bottom in step 6 is b Update based on the following formula: Where N is the effective pressure of the glacier, C is the maximum ratio of the shear stress at the bottom of the glacier to the effective pressure, and A is the maximum value of the ratio of the shear stress at the bottom of the glacier to the effective pressure. s is the glacier sliding parameter, B is the mobility parameter, n is the Glen flow law exponent, λ is the wavelength of the ice bed topography; γ is the roughness of the ice bed.

6. The ice thickness inversion method based on improved laminar flow of glacier bottom sliding according to claim 5 is characterized in that: In step 6, the effective pressure N of the glacier is based on the following formula: N=ρgH·cosα Where ρ represents the density of the ice body, g is the acceleration of gravity, H is the thickness of the glacier, and α is the slope of the glacier surface; The maximum value of the ratio of shear stress to effective pressure at the bottom of the glacier, C, is based on the following formula: C=0.84±0.02m max Among them, m max To obtain the maximum slope of the subglacial terrain based on the initial subglacial terrain Z0; The glacier sliding parameter is based on the following formula: Wherein, B is the fluidity parameter, B = 430 ± 40 MPa -3 yr -1 , λ is the wavelength of the ice bed topography, r is the roughness of the ice bed; The roughness of the ice bed is based on the following formula: Among them, k represents the number of pixels of the ice bed, j represents the sequence number of the pixel, and z j represents the elevation of the jth pixel of the ice bed, Indicates the average elevation of the ice bed.

7. A computer device comprising a memory and a processor, wherein the memory stores a computer program, wherein: When the processor executes the computer program, the steps of the ice thickness inversion method according to any one of claims 1 to 6 are implemented.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the ice thickness inversion method according to any one of claims 1 to 6 are implemented.

9. A computer program product comprising a computer program, characterized in that When the computer program is executed by a processor, the steps of the ice thickness inversion method according to any one of claims 1 to 6 are implemented.

Images