Ice thickness inversion method of improved laminar flow based on glacier bottom sliding

Abstract

The invention discloses a glacier bottom sliding-based improved laminar flow ice thickness inversion method, which comprises the following steps of: preparing grid-format glacier surface flow velocity, a digital elevation model (DEM) and vector-format glacier boundary data, calculating a glacier surface slope according to the digital elevation model (DEM), generating a glacier mask according to the glacier boundary data, and calculating the glacier surface slope according to the glacier surface slope. Unifying the five kinds of data into a projection coordinate system; constructing a laminar flow theoretical formula for calculating the glacier thickness H; assuming that the sliding speed of the glacier bottom is 0, and obtaining the initial glacier thickness; updating the valley shape factor, the bottom shear stress and the glacier bottom sliding speed in sequence to obtain a new glacier thickness; and calculating a difference value of adjacent glacier thicknesses, ending calculation and outputting a result when the difference value meets a set condition, otherwise, assigning the difference value to the initial glacier thickness, and restarting calculation until the difference value meets the set condition. According to the invention, the precision of the ice thickness model based on the traditional laminar flow theory is effectively improved, and the terrain under ice can be automatically obtained.

Description

Technical Field

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

[0002] Glaciers are natural ice bodies that exist in polar and alpine regions for a long time, formed by solid precipitation such as snow and hail through repeated processes such as compaction, crystallization, and freezing. The formation of mountain glaciers is jointly determined by climatic conditions such as precipitation and temperature, and the topographic conditions of the altitude of the mountains above the snow line, which directly affect the number and size of glaciers. The global glacier area is about 16 million square kilometers, accounting for about 11% of the total land area, and consolidating 68.7% of the world's freshwater resources. The world's mountain glaciers are mainly distributed in the Alps, Alaska, the high mountains of Asia, the southern Andes, and other regions, flowing slowly from high to low along the valley. As an important part of the cryosphere, glaciers are not only an indicator of climate change, but also one of the important factors restricting human socio-economic development. Due to the continuous global warming in recent decades, glaciers in alpine Asia have experienced severe recession and mass loss, leading to accelerated sea level rise. In addition, the rapid movement of glaciers may trigger various glacial disasters, such as glacial lake outburst floods (GLOFs) and ice avalanches. Therefore, developing ice thickness inversion models and clarifying the thickness distribution and volume of mountain glaciers are important foundations for glaciological research.

[0003] At present, regional and global scale mountain glacier thickness and volume estimation mainly uses the topographic characteristics, motion characteristics, and geometric characteristics of glaciers to invert through the principles of glacier dynamics. However, there are significant differences in the inversion of ice volume from existing ice thickness models based on laminar theory, which are attributed to the accuracy of ice thickness models, input data quality, and parameter estimation. In these models, laminar models estimate ice thickness based on glacier surface velocity and basal sliding. Generally, laminar models assume a scaling relationship between glacier basal sliding and surface velocity when estimating ice thickness. For example, it can be assumed that the ratio of basal sliding to surface velocity is 25%, that is, the ratio of slope to surface velocity, or that basal sliding is 0 in winter. However, basal sliding is affected by factors such as subglacial topography, basal shear stress, and ice cover. The ratio of basal sliding to surface velocity varies for different glaciers, and the measured ratios range from 0.03 to 1. Therefore, considering glacier basal sliding in laminar theory is crucial to improve the accuracy of ice thickness estimation. Glacier surface velocity can be regarded as the sum of basal sliding and ice deformation. Glacier basal sliding is closely related to the basal temperature and water content of the ice-rock interface. Generally, the subglacial terrain is irregular and has many obstacles. Glacier movement is driven by two mechanisms: pressure melting phenomenon and creep enhancement. However, obtaining information about the geometry and temperature of the ice bed remains challenging. Due to the technical complexity of field measurements, the slip law is usually used to simulate basal sliding in glacier dynamics research. In 1957, Veltman established the sliding law, which assumes that the glacier rests on hard bedrock, there is no shear stress and cavities at the interface, and the ice bed is idealized as a cubic bulge. However, the basal shear stress, effective pressure (the difference between the water pressure of the overlying layer of the ice layer and the base) and the roughness of the ice bed all affect the basal sliding. Therefore, the bed cavities need to be considered in the sliding law, and constraints are added to prevent the ratio of basal shear stress to effective stress from increasing with the increase of sliding speed. The Coulomb friction law satisfies the above conditions, so that the ratio of basal shear stress to effective stress first increases to the maximum value, then gradually decreases and tends to be stable. Therefore, the sliding law is used to determine the basal sliding of the laminar model to improve the accuracy of ice thickness estimation, but further research is lacking. Summary of the invention

[0004] The purpose of the present invention is to address the problem that the existing ice thickness model based on laminar flow theory does not consider the bottom sliding, resulting in low ice thickness estimation accuracy, and propose an improved laminar flow inversion ice thickness method based on glacier bottom sliding, which is achieved by combining the glacier bottom sliding law. The present invention can be used to estimate the bottom sliding of the glacier in space, and then calculate the glacier thickness through laminar flow theory, which helps to improve the ice thickness inversion accuracy, and can also automatically obtain the subglacial terrain, providing an important data basis 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 bottom sliding of glaciers 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 obtain the initial glacier thickness H using the glacier thickness calculation formula constructed in step 2 0 , and the initial glacier thickness H is subtracted from the digital elevation model DEM 0 Get the initial subglacial terrain Z 0 ;

[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 The 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 between two adjacent calculations of glacier thickness ΔH = H i -H i-1 , i≥1, i is the calculation number of glacier thickness. 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 H 0 , 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, the glacier boundary data, the glacier surface slope α, and the 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 They represent the glacier surface velocity and the glacier bottom sliding velocity respectively, n represents the Glen flow law index; A represents the ice creep coefficient; τ b represents the shear stress at the bottom of the glacier;

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

[0020] τ b =fρgHsinα

[0021] Among them, f represents the valley shape factor, ρ represents the density of the ice body, g is the gravitational acceleration, 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, 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 calculated is taken to obtain the average width w of the glacier.

[0029] The sliding speed u of the bottom of the glacier 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 value of the ratio of the shear stress at the bottom of the glacier to the effective pressure, and A is s is the glacier sliding parameter, B is the fluidity 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] Among them, ρ 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 is the initial subglacial terrain Z 0 Obtain the maximum slope of the subglacial terrain;

[0038] The glacier sliding parameters are 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, Represents the average elevation of the ice bed.

[0044] A computer device comprises a memory and a processor, wherein the memory stores a computer program, and wherein the processor implements the steps of the ice thickness inversion method described in any one of claims 1 to 7 when executing the computer program.

[0045] 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 described in any one of claims 1 to 7 are implemented. 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 described in any one of claims 1 to 7 are implemented.

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

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

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

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

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

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

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

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

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

[0055] 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 implementation examples described herein are only used to illustrate and explain the present invention and are not used to limit the present invention.

[0056] Embodiment 1:

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

[0058] Step 1: Obtain the glacier surface velocity u in raster format s and digital elevation model (DEM), obtain glacier boundary data in vector format. Calculate the glacier surface slope α based on the digital elevation model DEM; generate a glacier mask based on the glacier boundary data. s , digital elevation model, glacier boundary data, glacier surface slope, and glacier mask are unified into a projection 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.

[0059] 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 .

[0060] 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:

[0061]

[0062] Where: H (m) represents the thickness of the glacier; u s and u b are the surface velocity of the glacier and the sliding velocity at the bottom 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 temperature, structure and water content of the ice, and is 2.4×10 -24 Pa -3 s -1 ; τ b represents the shear stress at the bottom of the glacier, calculated using the shallow ice approximation formula.

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

[0064] 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 calculation formula for the glacier thickness H.

[0065]

[0066] Step 3: First, slide at the bottom of the glacier at a speed u b In the case of unknown, assume that the sliding speed u at the bottom of the glacier b The initial glacier thickness H is obtained by the glacier thickness calculation formula constructed in step 2. 0 , and the initial glacier thickness H is subtracted from the digital elevation model DEM 0 Get the initial subglacial terrain Z 0 .

[0067] 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.

[0068]

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

[0070] The average width w of the glacier is obtained as follows: extract the points on the center line of the glacier according to the glacier mask and calculate the coordinates of each point; then calculate twice the distance from each point to the glacier boundary in the direction perpendicular to the center line of the glacier as the glacier width, and take the average value to obtain the average width w of the glacier.

[0071] 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).

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

[0073]

[0074] Where N is 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 the ratio of shear stress to effective pressure at the bottom of the glacier can reach, C = 0.84 ± 0.02m max , where m max represents the maximum slope of the subglacial terrain, according to the initial subglacial terrain Z 0 The maximum slope m of the subglacial terrain is obtained 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), In the formula, 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, 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 .

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

[0076] 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 substituted into the glacier thickness calculation formula in step 2 to obtain the new glacier thickness H i .

[0077] Step 8: Calculate the difference between two adjacent calculations of glacier thickness Δ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 be assigned to the new initial glacier thickness H 0 , and the sliding speed u of the glacier bottom updated according to step (6) b , return to step 3.

[0078] Figure 5 The thickness distribution of ChhotaShigri glacier calculated by this method, Figure 6 Final subglacial topography of the ChhotaShigri Glacier calculated for this method.

[0079] A person of ordinary skill in the art can 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.

[0080] Embodiment 2:

[0081] 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-mentioned embodiment 1 when executing the computer program.

[0082] Embodiment 3:

[0083] 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.

[0084] Embodiment 4:

[0085] 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.

[0086] It should be noted that the embodiments described in the present invention are merely examples of the spirit of the present invention. Those skilled in the art may make various modifications or additions to the described embodiments or replace them in similar ways, but they will not deviate from the spirit of the present invention or exceed the scope defined by the attached claims.

Claims

1. An improved laminar flow inversion method for ice thickness based on bottom sliding of glaciers, characterized in that: 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 The value is assigned as 0m / s, and the initial glacier thickness H0 is obtained by the glacier thickness calculation formula 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; 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 The 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 between two consecutive calculations of glacier thickness ΔH = H i -H i-1 , i≥1, i is the calculation number of glacier thickness. 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 as the final subglacial terrain. i Return to step 3 using H0 as the new initial 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, the glacier boundary data, the glacier surface slope α, and the 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: 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 They represent the glacier surface velocity and the glacier bottom sliding velocity respectively, n represents the Glen flow law index; A represents the ice creep coefficient; τ b represents the shear stress at the bottom of the glacier; Glacier bottom shear stress τ b Based on the following formula: t b =fρgHsina Among them, f represents the valley shape factor, ρ represents the density of the ice body, g is the gravitational acceleration, 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:

4. The ice thickness inversion method based on improved laminar flow of glacier bottom sliding according to claim 3 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.

5. The ice thickness inversion method based on improved laminar flow of glacier bottom sliding according to claim 4 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, 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 calculated is taken to obtain the average width w of the glacier.

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

7. The ice thickness inversion method based on improved laminar flow of glacier bottom sliding according to claim 6 is characterized in that: In step 6, the effective pressure N of the glacier is based on the following formula: N=ρgH·cosα Among them, ρ 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 according to the initial subglacial terrain Z0; The glacier sliding parameters are 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, Represents the average elevation of the ice bed.

8. 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 7 are implemented.

9. 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 7 are implemented.

10. 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 7 are implemented.

Images