Coverage for C:\src\imod-python\imod\prepare\topsystem\resistance.py: 15%

27 statements  

« prev     ^ index     » next       coverage.py v7.5.1, created at 2024-05-08 14:15 +0200

1import numpy as np 

2 

3from imod.typing import GridDataArray 

4 

5 

6def c_radial( 

7 L: GridDataArray, 

8 kh: GridDataArray, 

9 kv: GridDataArray, 

10 B: GridDataArray, 

11 D: GridDataArray, 

12) -> GridDataArray: 

13 """ 

14 Ernst's radial resistance term to a drain. 

15 

16 Parameters 

17 ---------- 

18 L : xr.DataArray of floats 

19 distance between water features 

20 kh : xr.DataArray of floats 

21 horizontal conductivity 

22 kv : xr.DataArray of floats 

23 vertical conductivity 

24 B : xr.DataArray of floats 

25 water feature wetted perimeter 

26 D : xr.DataArray of floats 

27 saturated thickness of the top system 

28 

29 Returns 

30 ------- 

31 radial_c : xr.DataArray 

32 Ernst's radial resistance for a drain 

33 """ 

34 # Take anisotropy into account fully 

35 c = ( 

36 L 

37 / (np.pi * np.sqrt(kh * kv)) 

38 * np.log((4.0 * D) / (np.pi * B) * np.sqrt(kh / kv)) 

39 ) 

40 c = c.where(~(c < 0.0), other=0.0) 

41 return c 

42 

43 

44def c_leakage( 

45 kh: GridDataArray, 

46 kv: GridDataArray, 

47 D: GridDataArray, 

48 c0: GridDataArray, 

49 c1: GridDataArray, 

50 B: GridDataArray, 

51 length: GridDataArray, 

52 dx: GridDataArray, 

53 dy: GridDataArray, 

54) -> GridDataArray: 

55 """ 

56 Computes the phreatic leakage resistance. 

57 

58 Parameters 

59 ---------- 

60 kh : xr.DataArray of floats 

61 horizontal conductivity of phreatic aquifer 

62 kv : xr.DataArray of floats 

63 vertical conductivity of phreatic aquifer 

64 c0 : xr.DataArray of floats 

65 hydraulic bed resistance of water feature 

66 c1 : xr.DataArray of floats 

67 hydraulic resistance of the first aquitard 

68 D : xr.DataArray of floats 

69 saturated thickness of the top system 

70 B : xr.DataArray of floats 

71 water feature wetted perimeter 

72 length : xr.DataArray of floats 

73 water feature length per cell 

74 dx : xr.DataArray of floats 

75 cellsize in x 

76 dy : xr.DataArray of floats 

77 cellsize in y 

78 

79 Returns 

80 ------- 

81 c_leakage: xr.DataArray of floats 

82 Hydraulic resistance of water features corrected for intra-cell 

83 hydraulic resistance and surface water interaction. 

84 """ 

85 

86 def f(x): 

87 """ 

88 two x times cotangens hyperbolicus of x 

89 """ 

90 # Avoid overflow for large x values 

91 # 20 is safe for 32 bit floats 

92 x = x.where(~(x > 20.0), other=20.0) 

93 exp2x = np.exp(2.0 * x) 

94 return x * (exp2x + 1) / (exp2x - 1.0) 

95 

96 # TODO: raise ValueError when cells are far from square? 

97 # Since method of average ditch distance will not work well 

98 # Compute geometry of ditches within cells 

99 cell_area = abs(dy * dx) 

100 wetted_area = length * B 

101 fraction_wetted = wetted_area / cell_area 

102 # Compute average distance between ditches 

103 L = (cell_area - wetted_area) / length 

104 

105 # Compute total resistance to aquifer 

106 c1_prime = c1 + (D / kv) 

107 # Compute influence for the part below the ditch 

108 labda_B = np.sqrt((kh * D * c1_prime * c0) / (c1_prime + c0)) 

109 # ... and the field 

110 labda_L = np.sqrt(c1_prime * kh * D) 

111 

112 x_B = B / (2.0 * labda_B) 

113 x_L = L / (2.0 * labda_L) 

114 

115 # Compute feeding resistance 

116 c_rad = c_radial(L, kv, kh, B, D) 

117 c_L = (c0 + c1_prime) * f(x_L) + (c0 * L / B) * f(x_B) 

118 c_B = (c1_prime + c0) / (c_L - c0 * L / B) * c_L 

119 # total feeding resistance equals the harmonic mean of resistances below 

120 # ditch (B) and field (L) plus the radical resistance 

121 # Can also be regarded as area weighted sum of conductances of B and L 

122 c_total = 1.0 / (fraction_wetted / c_B + (1.0 - fraction_wetted) / c_L) + c_rad 

123 # Subtract aquifer and aquitard resistance from feeding resistance 

124 c = c_total - c1_prime 

125 # Replace areas where cells are fully covered by water 

126 c = c.where(~(L == 0.0), other=c0) 

127 return c