Code for "Multilayer Analytic Element Modeling of Radial Collector Wells"

[bibek-ux]

Hello team,
I have been exploring this repo for some research purposes. I studied the paper titled "Multilayer Analytic Element Modeling of Radial Collector Wells". The paper was very insightful and I have been trying to recreate the results from the paper. However I did not find any implementation of friction factor in the linesink elements code associated with the paper are not uploaded here. Can I go through those codes if it isn't confidential. I assure you if our work makes its way to publication, the respective authors or team will be thoroughly acknowledged and cited.

[mbakker7]

I am sorry to report that this code is not longer available as it was based on a very old version of TimML (well before TimML was even merged to github). Thank you for your interest.

[bibek-ux]

Hello team,
I went through the paper "Multilayer Analytic Element Modeling of Radial
Collector Wells" and was trying to replicate the results. The following is the implementation I came up with considering no entry loss or skin friction. However, keep on getting overestimation of discharge of about 31000m^3/d, ( 22400 in the paper). Any suggestions here ?

import numpy as np
import matplotlib.pyplot as plt
import timml as tml
from dataclasses import dataclass
from typing import List, Optional

Global parameters
G = 9.81 # gravity m/s²

https://github.com/DataClass
class ArmSegment:
"""Data structure for a single collector well arm segment."""
x1: float
y1: float
x2: float
y2: float
xc: float # center x
yc: float # center y
length: float

--- mutable state ---

strength: float
target_head: float = 20.0 # target head inside segment (m)
head_loss : float = 0.2

--- timml object ---

timml_element: Optional[tml.HeadLineSink] = None

class CollectorWellArm:
"""
A single arm of a collector well.
Manages its segments and calculates internal head losses.
"""
def init(self, x1: float, y1: float, x2: float, y2: float,
n_segments: int = 10, radius: float = 0.15, friction_factor: float = 0.02):
self.x1, self.y1 = x1, y1
self.x2, self.y2 = x2, y2
self.n_segments = int(n_segments)
self.radius = radius
self.f = friction_factor
self.length = np.sqrt((x2 - x1)**2 + (y2 - y1)**2)

Segments are ordered from tip (index 0) to caisson (index n-1)

self.segments = self._create_segments()

def _create_segments(self) -> List[ArmSegment]:
segments = []
seg_length = self.length / self.n_segments
x_p = np.linspace(self.x1, self.x2, self.n_segments + 1)
y_p = np.linspace(self.y1, self.y2, self.n_segments + 1)

  for i in range(self.n_segments):
      segments.append(ArmSegment(
          x1=float(x_p[i]), y1=float(y_p[i]), x2=float(x_p[i+1]), y2=float(y_p[i+1]),
          xc=float((x_p[i] + x_p[i+1]) / 2.0),
          yc=float((y_p[i] + y_p[i+1]) / 2.0),
          length=float(seg_length),
          strength = 100.0,  # initial guess (m^2/d)
          target_head = None,
          head_loss = 0.0
      ))
  # Set an initial target_head (will be updated later)
  for seg in segments:
      seg.target_head = 0.0
  return segments

def _get_cumulative_inflow(self, segment_index: int) -> float:
total_inflow = 0.0
# sum from tip (0) up to and including segment_index
for i in range(segment_index + 1):
total_inflow += self.segments[i].strength * self.segments[i].length
return total_inflow

def calculate_head_losses(self, caisson_head: float):
"""
Updates target_head for each segment using friction and exit loss,
working backwards from the caisson.
"""
m_last = self.n_segments - 1
# initialize last segment head (center of last segment)
exit_loss = self._calculate_exit_loss()
# head inside last segment center (close to caisson)
head_last = caisson_head + exit_loss
self.segments[m_last].head_loss = exit_loss
self.segments[m_last].target_head = head_last

  # Work backwards from the segment just upstream of last
  for m in range(m_last - 1, -1, -1):
      friction = self._calculate_friction_loss(m)
      head_last += friction
      self.segments[m].head_loss = friction
      self.segments[m].target_head = head_last

def calculate_friction_loss(self, m: int) -> float:
"""
Eq. 12 based friction loss between centers of segment m and m+1.
"""
# compute average discharge Qm (m^2/d length*strength) then convert to m^3/s
# Sm = sum{j=0}^{m-1} r_j L_j + (r_m L_m)/2 (discharge at center m)
# Qm = (Sm + Sm_prev)/2
# Qm = sum_{j=0}^{m-1} r_j L_j + (3r_mL_m + r_{m+1}*L_{m+1})/4

  if m < 0 or m >= self.n_segments - 1:
      return 0.0

  sum_part = 0.0
  for j in range(0, m):
      sum_part += self.segments[j].strength * self.segments[j].length

  rm = self.segments[m].strength
  rmp1 = self.segments[m+1].strength
  Lm = self.segments[m].length
  Lmp1 = self.segments[m+1].length

  Qm_m3_per_day = sum_part + (3.0 * rm * Lm + rmp1 * Lmp1) / 4.0  # m^3/d
  if Qm_m3_per_day <= 0.0:
      return 0.0

  Qm_s = Qm_m3_per_day / 86400.0  # m^3/s
  lm = (Lm + Lmp1) / 2.0  # distance between centers approximated as average segment length
  # Darcy-Weisbach friction (paper eq.10)
  delta_phi = (self.f * lm * Qm_s**2) / (4.0 * G * np.pi**2 * (self.radius)**5)
  return float(delta_phi)

def _calculate_exit_loss(self) -> float:
"""
Exit loss at the last segment (paper eq.15) using QM from eq.16.
"""
m_last = self.n_segments - 1
if m_last < 0:
return 0.0

  sum_part = 0.0
  for j in range(0, m_last):
      sum_part += self.segments[j].strength * self.segments[j].length

  rM = self.segments[m_last].strength
  LM = self.segments[m_last].length
  QM_m3_per_day = sum_part + (3.0 * rM * LM) / 4.0
  if QM_m3_per_day <= 0.0:
      return 0.0

  QM_s = QM_m3_per_day / 86400.0
  lM = LM / 2.0  # distance between center and caisson approx half last segment length
  delta_phi_M = (self.f * (lM / (2.0 * self.radius)) + 1.0) * \
                (QM_s**2 / (2.0 * G * (np.pi * self.radius**2)**2))
  return float(delta_phi_M)

def get_total_arm_discharge(self) -> float:
"""Get total inflow into this arm by summing all segment inflows (Eq. 7)."""
return self._get_cumulative_inflow(self.n_segments - 1)
class CollectorWell:
"""
Collector Well class that manages all arms and the iterative TimML solution.
"""
def init(self, model: tml.Model, xc: float = 0, yc: float = 0, caisson_radius: float = 3,
n_arms: int = 3, arm_length: float = 60, arm_radius: float = 0.15,
friction_factor: float = 0.02, n_segments_per_arm: int = 10,
layer: int = 13, caisson_head: float = 20):
self.model = model
self.xc, self.yc = xc, yc
self.caisson_radius = caisson_radius
self.caisson_head = caisson_head
self.arms = self._create_arms(n_arms, arm_length, arm_radius,
friction_factor, n_segments_per_arm, layer)
self.layer = layer

def _create_arms(self, n_arms, arm_length, arm_radius, friction_factor, n_segments, layer) -> List[CollectorWellArm]:
arms = []
angles = np.linspace(0, 2 * np.pi, n_arms, endpoint=False)
for angle in angles:
x1 = self.xc + self.caisson_radius * np.cos(angle)
y1 = self.yc + self.caisson_radius * np.sin(angle)
x2 = self.xc + (self.caisson_radius + arm_length) * np.cos(angle)
y2 = self.yc + (self.caisson_radius + arm_length) * np.sin(angle)
arms.append(CollectorWellArm(x1, y1, x2, y2, n_segments, radius=arm_radius, friction_factor=friction_factor))
return arms

def solve(self, max_iterations: int = 50, tolerance: float = 1e-2):
"""
Solves the model iteratively to account for non-linear friction losses.
"""
print("Starting iterative solution for collector well...")
print(f"Target caisson head: {self.caisson_head} m")
print("=" * 60)

  for i in range(max_iterations):
      # 1. Setup/Update TimML model with HeadLineSinks
      self._setup_timml_model()

      # 2. Solve the linear groundwater flow problem
      self.model.solve()

      # 3. Update strengths from the TimML solution
      self._update_strengths_from_timml()

      # 4. Calculate new target heads for the next iteration
      self._update_target_heads()
      
      # 5. Check for convergence
      total_discharge = self.get_total_discharge()
      if i > 0:
          change = np.abs(total_discharge - old_total_discharge)
          print(f"Iteration {i+1:2d}: Total Q = {total_discharge:7.1f} m³/d, Change = {change:7.2f} m³/d")
          if change < tolerance:
              print(f"\nConverged after {i+1} iterations.")
              break
      else:
          print(f"Iteration {i+1:2d}: Total Q = {total_discharge:7.1f} m³/d")

      old_total_discharge = total_discharge
  else:
      print("\nWarning: Maximum iterations reached without convergence.")
  
  print(f"Final total discharge: {self.get_total_discharge():.1f} m³/d")

def _setup_timml_model(self):
# Remove old elements
for arm in self.arms:
for seg in arm.segments:
if seg.timml_element is not None and seg.timml_element in self.model.elementlist:
self.model.remove_element(seg.timml_element)
seg.timml_element = None

  # Add new elements
  for arm_idx, arm in enumerate(self.arms):
      for seg_idx, seg in enumerate(arm.segments):
          # Determine the head at the downstream end of this segment:
          if seg_idx < arm.n_segments - 1:
              head_end = arm.segments[seg_idx + 1].target_head
          else:
              # last segment -> downstream is caisson head
              head_end = self.caisson_head

          head_start = seg.target_head if seg.target_head is not None else self.caisson_head

          ls = tml.HeadLineSinkString(self.model,
                                xy = [(seg.x1, seg.y1), (seg.x2, seg.y2)],
                                hls=[head_start, head_end],
                                layers=[self.layer],  # ensure list
                                res = 0,
                                wh = 0.00005,
                                label=f'Arm{arm_idx}_Seg{seg_idx}',
                                order = 3,
                               )
          seg.timml_element = ls

def _update_strengths_from_timml(self, relaxation=0.5):
for arm in self.arms:
for seg in arm.segments:
discharge = seg.timml_element.discharge()[self.layer]

          strength = seg.timml_element.parameters[:, 0]
          # print(f"Segment strength: {strength} m²/d, Discharge: {discharge:.2f} m³/d")
         

          seg.strength = seg.strength * (1 - relaxation) + relaxation * strength[0]

def _update_target_heads(self):
"""Update target heads for all segments based on current friction losses."""
for arm in self.arms:
print(f"Target head: {arm.segments[0].target_head:.2f} m, Head loss: {arm.segments[0].head_loss:.6f} m")
arm.calculate_head_losses(self.caisson_head)

def get_total_discharge(self) -> float:
"""Calculate the total discharge of the entire collector well."""
return sum(arm.get_total_arm_discharge() for arm in self.arms)

def plot_results(self):
"""Generates four plots to visualize the converged solution."""
print("Generating result plots...")
plt.figure(figsize=(14, 10))
plt.suptitle('Collector Well Performance Analysis', fontsize=16)

  # Plot 1: Discharge distribution (strength)
  ax1 = plt.subplot(2, 2, 1)
  for i, arm in enumerate(self.arms):
      distances = [(arm.n_segments - j - 0.5) * arm.segments[0].length for j in range(arm.n_segments)]
      strengths = [seg.strength for seg in arm.segments]
      ax1.plot(distances, strengths, 'o-', label=f'Arm {i+1}')
  ax1.set_xlabel('Distance from Caisson (m)')
  ax1.set_ylabel('Inflow Strength (m²/d)')
  ax1.set_title('Discharge Distribution Along Arms')
  ax1.legend()
  ax1.grid(True, alpha=0.5)

  # Plot 2: Head distribution
  ax2 = plt.subplot(2, 2, 2)
  for i, arm in enumerate(self.arms):
      distances = [(arm.n_segments - j - 0.5) * arm.segments[0].length for j in range(arm.n_segments)]
      heads_inside = [seg.target_head for seg in arm.segments]
      print(heads_inside)
      ax2.plot(distances, heads_inside, 'o-', label=f'Arm {i+1}')
  ax2.axhline(self.caisson_head, color='k', ls='--', label='Caisson Head')
  ax2.set_xlabel('Distance from Caisson (m)')
  ax2.set_ylabel('Head Inside Arm (m)')
  ax2.set_title('Head Distribution Inside Arms')
  ax2.legend()
  ax2.grid(True, alpha=0.5)

  # Plot 3: Cumulative discharge
  ax3 = plt.subplot(2, 2, 3)
  for i, arm in enumerate(self.arms):
      distances = [(arm.n_segments - j - 0.5) * arm.segments[0].length for j in range(arm.n_segments)]
      cum_q = [arm._get_cumulative_inflow(j) for j in range(arm.n_segments)]
      ax3.plot(distances, cum_q, 'o-', label=f'Arm {i+1}')
  ax3.set_xlabel('Distance from Caisson (m)')
  ax3.set_ylabel('Cumulative Discharge (m³/d)')
  ax3.set_title('Cumulative Discharge Along Arms (Tip to Caisson)')
  ax3.legend()
  ax3.grid(True, alpha=0.5)
  
  # Plot 4: Aquifer head contour
  ax4 = plt.subplot(2, 2, 4)
  plt.sca(ax4)
  self.model.contour(win=[-30, 90, -60, 60], ngr=60, layers=self.layer,
                     levels=15, labels=True , decimals=3)
  ax4.set_aspect('equal')
  ax4.set_title(f'Head Contours in well Layer  (m)')
  
  plt.tight_layout(rect=[0, 0.03, 1, 0.95])
  plt.show()

def run_simulation():
"""Sets up and runs the collector well simulation."""
layer_tops = [24, 16, 12, 10, 8, 7, 6, 5.45, 5.15, 4.85, 4.55, 4, 3.45, 3.15, 2.85, 2.55, 1.95, 1, 0]
kaq = np.append(100 * np.ones(3), 200 * np.ones(15))
kzoverkh = np.append((60/100) * np.ones(3), (60/200) * np.ones(15))
layer_idx = 13

  # Create TimML model

ml = tml.Model3D(kaq=kaq, z=layer_tops,kzoverkh = kzoverkh , topboundary='conf', hstar=0)

tml.Constant(ml, xr=200, yr=0, hr=24 , layer = 0)

cw = CollectorWell(
model=ml,
n_arms=3,
arm_length=60,
n_segments_per_arm=10,
caisson_head=20,
friction_factor=0.02,
layer=layer_idx, # will compute layer_idx below to match arm elevation
arm_radius=0.15, # match paper (diameter = 0.3 m)
caisson_radius=3
)

Solve and visualize

cw.solve(max_iterations=200, tolerance=1e-4)
print(f"head at top is {ml.head(0,0)[0]}")
print(f"Head at end point in collecter arm is {ml.head(60 , 0)[layer_idx]}")
print(f"Head at caisson end point in collecter arm is {ml.head(6 , 0)[layer_idx]}")

cw.plot_results()
if name == "main":
run_simulation()

[mbakker7]

I think the main reason you were not able to reproduce the multi-layer results of the Bakker, Kelson, Luther (2005) paper is that in TimML the control points are placed on the line-sinks, while in the paper the control points are placed away from the line-sinks at a distance equal to the radius of the arms.

I made a few additions on the controlpoints_off_linesink branch to allow the user to specify the distance between the control points and the line-sink (the dely keyword). I also modified the line-sinks string elements to allow for the specification of a non-continuous string (as is the case for a collector well with multiple arms) by specifying xy as an array with 4 columns, where each row is x1, y1, x2, y2. Finally, I added a function headinside to the line-sinks that can be used to return the head inside the line-sink. Mind you: this last function does not yet work for line-sinks with an entry resistance. I created an issue #135 to add that later (feel free to submit a PR if you want to).

In the notebook bakker_kelson_luther_2005 click here I simulated the first two examples of the paper and get the same results. For the case with flow resistance in the arms, you have to iterate (as mentioned in the paper and in the script you submitted above). Feel free to add iteration to the above notebook and submit a PR.

Thanks for bringing this to our attention. I hope this helps.