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

26 statements  

« prev     ^ index     » next       coverage.py v7.4.4, created at 2024-04-08 10:26 +0200

1import numpy as np 

2 

3 

4def c_radial(L, kh, kv, B, D): 

5 """ 

6 Ernst's radial resistance term to a drain. 

7 

8 Parameters 

9 ---------- 

10 L : xr.DataArray of floats 

11 distance between water features 

12 kh : xr.DataArray of floats 

13 horizontal conductivity 

14 kv : xr.DataArray of floats 

15 vertical conductivity 

16 B : xr.DataArray of floats 

17 water feature wetted perimeter 

18 D : xr.DataArray of floats 

19 saturated thickness of the top system 

20 

21 Returns 

22 ------- 

23 radial_c : xr.DataArray 

24 Ernst's radial resistance for a drain 

25 """ 

26 # Take anisotropy into account fully 

27 c = ( 

28 L 

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

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

31 ) 

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

33 return c 

34 

35 

36def c_leakage(kh, kv, D, c0, c1, B, length, dx, dy): 

37 """ 

38 Computes the phreatic leakage resistance. 

39 

40 Parameters 

41 ---------- 

42 kh : xr.DataArray of floats 

43 horizontal conductivity of phreatic aquifer 

44 kv : xr.DataArray of floats 

45 vertical conductivity of phreatic aquifer 

46 c0 : xr.DataArray of floats 

47 hydraulic bed resistance of water feature 

48 c1 : xr.DataArray of floats 

49 hydraulic resistance of the first aquitard 

50 D : xr.DataArray of floats 

51 saturated thickness of the top system 

52 B : xr.DataArray of floats 

53 water feature wetted perimeter 

54 length : xr.DataArray of floats 

55 water feature length per cell 

56 dx : xr.DataArray of floats 

57 cellsize in x 

58 dy : xr.DataArray of floats 

59 cellsize in y 

60 

61 Returns 

62 ------- 

63 c_leakage: xr.DataArray of floats 

64 Hydraulic resistance of water features corrected for intra-cell 

65 hydraulic resistance and surface water interaction. 

66 """ 

67 

68 def f(x): 

69 """ 

70 two x times cotangens hyperbolicus of x 

71 """ 

72 # Avoid overflow for large x values 

73 # 20 is safe for 32 bit floats 

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

75 exp2x = np.exp(2.0 * x) 

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

77 

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

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

80 # Compute geometry of ditches within cells 

81 cell_area = abs(dy * dx) 

82 wetted_area = length * B 

83 fraction_wetted = wetted_area / cell_area 

84 # Compute average distance between ditches 

85 L = (cell_area - wetted_area) / length 

86 

87 # Compute total resistance to aquifer 

88 c1_prime = c1 + (D / kv) 

89 # Compute influence for the part below the ditch 

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

91 # ... and the field 

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

93 

94 x_B = B / (2.0 * labda_B) 

95 x_L = L / (2.0 * labda_L) 

96 

97 # Compute feeding resistance 

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

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

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

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

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

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

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

105 # Subtract aquifer and aquitard resistance from feeding resistance 

106 c = c_total - c1_prime 

107 # Replace areas where cells are fully covered by water 

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

109 return c