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
« prev ^ index » next coverage.py v7.5.1, created at 2024-05-08 14:15 +0200
1import numpy as np
3from imod.typing import GridDataArray
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.
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
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
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.
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
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 """
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)
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
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)
112 x_B = B / (2.0 * labda_B)
113 x_L = L / (2.0 * labda_L)
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