Eine Grundwasserabsenkung soll bei Tiefbauarbeiten gewährleisten, dass in einer Baugrube trockenen Fußes gearbeitet werden kann. Dazu wird der Grundwasserstand unterhalb der Baugrubensohle abgesenkt(und noch sicherheitshalber etwas darüber hinaus). Man möchte einen stabilen Absenktrichter im Grundwasserleiter erzeugen. Wie man so eine Absenkung berechnet, soll hier gezeigt werden.
   
   
   


Zylinder
Wir stellen uns einen Zylinder vor, der bis zum Rand mit Wasser gefüllt ist. Radius und die Höhe des Zylinders seien 1:


#Parameter
radius=1
wasserstand_innen=1
Wir errechnen aus diesen beiden Werten das Volumen. Vorher importieren wir noch das "math"- Modul, damit wir die Zahl π nicht tippen müssen:

import math
#Parameter
radius=1
wasserstand_innen=1
grundflaeche = math.pi * (radius**2)

volumen=grundflaeche*wasserstand_innen


Lassen wir die berechneten Werte ausgeben, erhalten wir π für das Volumen. (π ist zwar eine irrationale Zahl, aber das Volumen ist endlich.)
Stellen wir uns vor, eine Pumpe sei in der Mark des inneren Zylinders, ganz nach Euklid als "Länge ohne Breite". Jetzt pumpen wir Wasser aus diesem Zylinder mit einer Rate von 0,5 Volumeneinheiten pro Zeiteinheit und schauen, wie viele Zeitschritte es braucht, bis er leer ist, also volumen= 0 wird. Dazu machen wir eine while-Schleife, die solange von unserem Anfangsvolumen 0.5 Volumeneinheiten subtrahiert wie volumen >0 ist. Also während (while) dieser oder jener Zustand währt: Etwas machen.

import math

#Parameter
radius=1
wasserstand_innen=1
grundflaeche = math.pi * (radius**2)
abpumprate= 0.5

volumen=grundflaeche*wasserstand_innen

while volumen>=0:
  print(f"Volumen: {volumen:.2f}")
  volumen -= abpumprate 

Beim 7. Zeitschritt wäre er also leer. Spielen wir das mal für verschiedene Radien durch, sehen wir, dass es mit zunehmendem Volumen bei gleicher Pumprate mehr Zeitschritte braucht, die entsprechenden Zylinder zu leeren. 

Stellen wir uns zwei weitere Zylinder mit gleichem Volumen als Ringe konzentrisch um den inneren Vollzylinder angeordnet vor. 
 
Berücksichtigen wir dabei auch, dass Ringe im Vergleich zu Vollzylindern bei gleichem Volumen einen größeren Radius benötigen. Wurzel(2) für den zweiten Zylinder und Wurzel(3). Die inneren Zylinderwände seien durchlöchert. 
Ein Mensch hat herausgefunden, dass das Volumen, das in einem Zeitschritt durch die Mantelfläche des inneren Zylinders geht, davon abhängig ist,
wie groß die Löcher sind,
wie groß die Differenz zwischen den Wasserständen in den Zylindern
und umgekehrt proportional von der Fließstrecke abhängig ist
. 
Als Formel notiert sieht dieser Zusammenhang so aus: Durchfluss = Mantelfläche*Lochweite*((Wasserstand_äußerer_zylinder-Wasserstand_innerer_zylinder)/Fließstrecke)
Die Mantelfläche der Zylinder kennen wir, den Durchfluss auch, denn er entspricht ja unserer Pumprate pro Zeitschritt. Die Lochweite sei 0,5, um diesen Wert ergänzen wir die berechneten Radien, damit das Volumen in allen Zylinder stets gleich bleibt. Die Fließstrecke durch die Wand sei 0.1 (Wanddicke)innerhalb eines Eimers ist das Wasser in einem freien Raum. Dort kann es sich ohne Widerstand bewegen, weswegen sie nicht zur Fließstrecke mitgezählt wird. Welche Wasserstände stellen sich an der Innenseite der Wände der inneren Zylinder ein? 

import math

abpumprate = 0.5
lochweite = 0.5
wanddicke = 0.1
anfangswasserstand = 1.0

# Wandradien fuer gleiche Volumina: r1 = √1, r2 = √2, r3 = √3
for i in range(1, 4):
    r = math.sqrt(i)

    mantelflaeche = 2 * math.pi * r * anfangswasserstand
    delta_h = (abpumprate * wanddicke) / (mantelflaeche * lochweite)
    wasserstand = anfangswasserstand - delta_h

    print(f"Wand {i} (r={r:.2f}) | Wasserstand innen: {wasserstand:.4f}")

Wir sehen zunächst, dass von Zylinder zu Zylinder der Wasserstand im ersten Zeitschritt von innen nach außen zunimmt. Bei ansonsten gleichen Parametern (Wanddicke, Abpumprate und Anfangswasserstand) sich mit kleinerer Lochweite. Die Mantelflächen werden immer größer, sodass der Wasserstand immer weniger abnehmen muss, um den gleichen Durchfluss zu erreichen. Lassen wir nun die Zylinder unendlich klein werden und entsprechend die Wanddicke unedlich ausdehnen, füllt das durchlöcherte Material den gesamten Raum lückenlos und wir erhalten Verhältnisse, wie sie im Boden herrschen. Den Löchern in unserem Beispielmodell ist dann z.B die Korngröße des Bodens analog. Je größer die Abpumprate, desto größer der Unterschied in den Wasserständen und umso größer der Durchfluss. Wenn wir diesmal nur die Lochweite variieren, sehen wir, dass kleinere Löcher zu einem geringeren Durchfluss führen. Da Löcher in einer Wand immer einen gewissen Strömungswiderstand besitzen, kann das Wasser aus dem äußeren Ring nicht beliebig schnell ins Zentrum nachfließen.

Der äußerste Zylinder habe nun einen konstanten Wasserstand, es ströme also ständig Wasser nach, sodass der Wasserstand stets dem Anfangswasserstand von 1 entspräche. Wie hoch müsste die Abpumprate sein, damit bei gegebener Durchlässigkeit der zweite Zylinder zur Hälfte gefüllt ist?



integriert über den Radius. was heißt das?

Umsetzung: Löse die Aufgabe wahlweise analytisch (durch Umstellen der Gleichung) oder numerisch (indem du dein Skript $Q$ so lange anpassen lässt, bis $h_2 \approx 0{,}5$ erreicht ist).










 

 Wellpoint-Verfahren
 Um an den Wänden einer z. B. $50 \times 50\text{ m}$ großen Baugrube die Ziel-Absenktiefe zu erreichen, musst du in der Mitte bei der Pumpe viel tiefer als nötig absenken. Du ziehst unzählige Kubikmeter Wasser ab, nur um die Ecken der Baugrube trocken zu bekommen.Ihre Trichter überlagern sich so, dass eine flache Wannenform entsteht, die exakt der Form der Baugrube entspricht. Es wird kein Wasser "zu tief" abgepumpt.
 
import math





 




import math
import matplotlib.pyplot as plt

# Parameter
radius_eimer1 = 1
radius_eimer2 = math.sqrt(5)
radius_eimer3 = math.sqrt(14)
höhe_eimer1 = 1
abpumprate = 0.5

# Grundflächen (A = pi * r^2 bzw. Ringflächen)
grundfläche_eimer1 = math.pi * (radius_eimer1**2)
grundfläche_eimer2 = math.pi * (radius_eimer2**2 - radius_eimer1**2)
grundfläche_eimer3 = math.pi * (radius_eimer3**2 - radius_eimer2**2)

# Startvolumina
volumen_eimer1 = grundfläche_eimer1 * höhe_eimer1
volumen_eimer2 = grundfläche_eimer2 * höhe_eimer1
volumen_eimer3 = grundfläche_eimer3 * höhe_eimer1

# Ausgewählte Zeitschritte zum Plotten
zeitschritte = [0, 1, 2, 3]

plt.figure(figsize=(9, 5))

for step in zeitschritte:
    # Verbleibendes Volumen nach 'step' Abpumpschritten
    v1 = max(0, volumen_eimer1 - step * abpumprate)
    v2 = max(0, volumen_eimer2 - step * abpumprate)
    v3 = max(0, volumen_eimer3 - step * abpumprate)
    
    # Wasserstände
    h1 = v1 / grundfläche_eimer1
    h2 = v2 / grundfläche_eimer2
    h3 = v3 / grundfläche_eimer3
    
plt.figure(figsize=(8, 5))

radii = [radius_eimer1, radius_eimer2, radius_eimer3]

for step in zeitschritte:
    v1 = max(0, volumen_eimer1 - step * abpumprate)
    v2 = max(0, volumen_eimer2 - step * abpumprate)
    v3 = max(0, volumen_eimer3 - step * abpumprate)
    
    h1 = v1 / grundfläche_eimer1
    h2 = v2 / grundfläche_eimer2
    h3 = v3 / grundfläche_eimer3
    
    plt.plot(radii, [h1, h2, h3], marker="o", label=f"Zeitschritt {step}")

plt.xlabel("Eimerradius [m]")
plt.ylabel("Wasserstand [m]")
plt.title("Wasserstand in Abhängigkeit vom Eimerradius")
plt.xticks(radii, labels=[f"r1={radius_eimer1:.2f}", f"r2={radius_eimer2:.2f}", f"r3={radius_eimer3:.2f}"])
plt.legend()
plt.grid(True)
plt.show()


Stellen wir uns vor, dass die Eimerwände durchlöchert seien.

 kf-Wert
Das Beispiel vollkommener Brunnen, der von oben bis unten hin verfiltert ist. 

Welche Abpumprate braucht es für eine bestimmte Absenkung in der Brunnenmitte? Welche Absenkung brauchen wir wenn (von oben betrachtet dann) eine Reichweite bis x m hergestellt werden soll.
Dann Superpositionsprinzip.
Dann unvollkommene Brunnen,also die Abnahme hängt auch vom Volumen des Brunnens ab und der Menge die pro Zeitrschritt gefördert werden kann. Außerdem Brunneneintrittsverluste, Anströmung von unten. Gespannte Grundwasserleiter.


 Die Wände dieser Eimer seien ganz fein durchlöchert, sodass die Abpumprate von 0.5 im ersten Eimer noch gilt und mit jeder Eimerwand 

Wenn wir Wasser fördern, muss es durch die Mantelfläche des Filterrohrs. Die Mantelfläche eines Zylinders errechnet sich wie folgt: A = 2*pi*r*h. Zeichnen wir um unseren Brunnenzylinder einen größeren Zylinder


   
   
     Oberer und unterer Querschnitt sind nicht durchflusswirksam und wir vernachlässigen einen Filterwiderstand.   Wir kennen den Ruhewasserpegel H und den gewünschten Wasserstand h. Je mehr wir absenken wollen, desto größer die Differenz der beiden, desto mehr Wasser muss also gefördert werden. Wir halten fest: Je länger die Filterstrecke, je größer die Absenkung, desto höher die nötige Förderrate. Je nachdem, ob unser Sandkasten aus Feinssand oder Grobsand besteht, wird die gewünschte Absenkung sich bei geringeren bzw. höheren Abflüssen einstellen. K gehört damit auch in unsere Formel. Alle diese Parameter gehen also direkt proportional mit einer höheren Förderrate einher, sprich dieser Zusammenhang wird durch einfache Multiplikation der Parameter ausgedrückt: Q=2*pi*h*K*(H-h). Diese Formel ist hinreichend genau für einen Sandkasten, der kaum größer ist als der Brunnen selbst.
  Je weiter man aber weggeht vom Brunnenmittelpunkt, desto geringer wird bei gegebener Querschnittsfläche und Leitfähigkeit der Abfluss. Weil wir ja aber mit gleicher Förderrate fördern, an dem Ruhewasserpegel nichts ändern können und immer noch das gleiche Absenkungsziel haben, die Durchlässigkeit des Bodens auch gleich bleibt und der Gradient mit zunehmender Enfernung kleiner wird, müssen immer größere Querschnitte durchströmt werden, damit das Kontinuitätsgesetz nicht verletzt wird.

Das heißt die Mantelfläche, die durchflossen wird, wird durch immer größere Zylinder um unseren Brunnen beschrieben. Wer in der Badewanne den Abfluss oder einen Gullideckel in den ein Hochwasser abfließt, beobachtet hat, weiß, dass sich ein Absenkungstrichter bildet, der mit zunehmender Entfernung zur Mitte hin immer flacher wird. Die Entfernung, bei der noch eine Absenkung stattfindet, wird Reichweite der Absenkung genannt. Es bildet sich ein Absenkungstrichter um den Brunnen. 

Durch jeden Fließquerschnitt (Mantelläche) muss diesselbe Wassermenge strömen. Das sagt das Kontinuitätsgesetz.
Rein prinzipiell: Welche Funktionsklasse beschreibt dieses Verhältnis am ehesten?

Radius-Ersatzbrunnen
Speicherkoeffizient

   
   
Brunnenformel nach Dupuit-Thiem
  
   Dupuit–Thiem method, which is an analytical model for the evaluation of steady-state pumping test
   
   To delineate the geometry of the cone of depression and the radius of influence, he installed several observation wells, both perpendic- ular and parallel to the estimated flow direction of ground- water
   
   
   In terms of dewatering, both fine and coarser sediments can be problematic. Fine and coarse sediments can be heavily saturated and difficult to dewater due to varying permeability or hydraulic conductivity. Generally speaking, finer sediments cannot sustain larger pumping rates required in dewatering projects. Coarse sediments require greater pumping rates and more demanding groundwater discharge management methods to adequately dewater higher hydraulic conductivity sediments. An understanding of the aquifer sediment size and sorting is crucial because increased groundwater velocities from aggressive pumping can mobilize finer sediments from within the aquifer and create piping (internal erosion) conditions that can lead to ground subsistence and/or pumping equipment destruction.
   
Welche Förderrate ist nötig, um mit einem Brunnen allein die gewünschte Absenkung zu erreichen?
Welche Förderraten ist nötig, um mit mehreren kleinen Brunnen die gewünschte Absenkung zu erreichen?
Wie viele Brunnen müssen dafür abgeteuft werden und wie sind diese zu positionieren?
Wann wird der gewünschte Absenkungsbetrag erreicht?
Wie groß ist die Reichweite der Grundwasserabsenkung? 
  
Werden Nutzungen oder Gebäude ggf. durch die Maßnahme beeinträchtigt?


   Die am häufigsten angewendeten Methoden zur Grundwasserabsenkung sind Förderbrunnen mit großem Durchmesser oder Spüllanzen, die um die Baugrube herum positioniert werden oder mehrere Brunnen kleineren Durchmessers (Punktbrunnenanlagen) im Wellpoint-Verfahren. Wird eine Vakuumpumpe genutzt, um Grundwasser zu fördern, beschränkt sich der maximal förderbare Höhenunterschied auf etwa 6m. Im Grunde drückt die Atmosphäre dabei das Wasser in die Schläuche, wo der Unterdruck herrscht. Die Atmosphäre hat einen Druck von 1 bar auf Meereshöhe. Das entspricht eigentlich 10,33 m Wassersäule, die ein Vakuum fördern können sollte; diese werden aber aufgrund von Reibungsverlusten, Kavitation und Luftausgasung in den Schläuchen nie erreicht.
    They
were spaced more closely to the well and further apart from it
   
   
    Vom Einfachen zum Schweren. Es gibt zwei Herangehensweisen. 
    Analytisch: das heißt eine mathematische Gleichung, die eine eindeutige Lösung für eine bestimmte Konstellation liefert.
    Numerisch: das heißt man teilt ein Gebiet in diskrete geometerische Objekte(Dreiecke, Quadrate, Rechtecke,...) und rechnet für jedes Objekt den Wasserstand aus, der sich einstellen müsste, wenn aus einer benachbarten Zelle so-und-so-viel ankommt und die abstromseitige Zelle so-und-so-viel ausfließt.
     Ich möchte hier beide Herangehensweisen für die Berechnung von Grundwasserabsenkungen nutzen und vergleichen.
      Ziel: Grundwasserabsenkung verstehen und berechnen so wie hier, hier, oder hier, aber mit Freeware wie LibreOffice Calc, MODFLOW oder Python

  
  Nachdem wir FloPy in Jupyter Notebook Startklar gemacht haben, erstellen wir unsere Baugrube. Zu allererst importieren wir alle nötigen Bibliotheken fürs Rechnen und Darstellen. Mit "os" können wir auf Ordner auf unseren Computern über die Pythonkonsole zugreifen. Mit numpy erstellen wir dataFrames. Das sind Daten im Tabellenformat. Mit pyplot aus dem matplotlib-Modul Visualisieren wir später unsere Ergebnisse. FloPy ist wie gesagt das MODFLOW-Ansprache Modul für Python.
  
import os
import numpy as np
import matplotlib.pyplot as plt
import flopy

Dann definieren wir den Speicherort für die MODFLOW-Dateien (ws, workspace), geben ihnen einen Namen und legen den Pfad mit dem MODFLOW-Kernel fest.

import os
import numpy as np
import matplotlib.pyplot as plt
import flopy

ws = "./baugrube_model"                            
sim_name = "baugrube"                             
exe_path = os.path.expanduser(r"PFAD/modflow_bin/mf6")


Jetzt kommen weitere Definitionen, nämlich über die Simulation(sim), Zeit-Diskreteisierung (TDIS), Gleichungs-Solver (IMS) und das Grundwassermodell.
TDIS erwartet einen Tupel mit folgenden Werten (Länge der Stressperiode "perlen, Anzahl der Time-Steps nstp, Multiplikator für Zeitlängen bei instationären Berechnungen  tsmult). In MODFLOW wird standardmäßig eine Periode ohne Angabe von Transienteneigenschaften als stationär (steady-state) interpretiert.
import os

import numpy as np
import matplotlib.pyplot as plt
import flopy

ws = "./baugrube_model"                            
sim_name = "baugrube"                             
exe_path = os.path.expanduser("~/modflow_bin/mf6")

sim = flopy.mf6.MFSimulation(sim_ws=ws, exe_name=exe_path)
tdis = flopy.mf6.ModflowTdis(sim, nper=1, perioddata=[(1.0, 1, 1.0)])
ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE")
gwf = flopy.mf6.ModflowGwf(sim, modelname=sim_name)



Als nächstes bauen wir die Geometrie, die unser Grundwasserleiter sein soll. Wir definieren dazu 
die Anzahl an Schichten (1), Zeilen (50) und Spalten (50)
die Dimensionen der Gitterzellen, also Zellbreite (delr) & Zelllänge (delc) in  (im Beispiel jeweils 10m)
die Geländehöhe (top) & Aquifersohle (btm)
die Durchlässigkeit (npf) in m/tag (z.B entsprechen 1e-4 m/s  8.64 m/d)
den Anfangswasserstand (ic) in m

import os

import numpy as np
import matplotlib.pyplot as plt
import flopy

ws = "./baugrube_model"                            
sim_name = "baugrube"                             
exe_path = os.path.expanduser(r"PFAD/modflow_bin/mf6")

sim = flopy.mf6.MFSimulation(sim_ws=ws, exe_name=exe_path)
tdis = flopy.mf6.ModflowTdis(sim, nper=1, perioddata=[(1.0, 1, 1.0)])
ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE")
gwf = flopy.mf6.ModflowGwf(sim, modelname=sim_name)

nlay, nrow, ncol = 1, 50, 50
dis = flopy.mf6.ModflowGwfdis(
    gwf,
    nlay=nlay, nrow=nrow, ncol=ncol,
    delr=10.0, delc=10.0,
    top=20.0, botm=0.0)

npf = flopy.mf6.ModflowGwfnpf(gwf, k=8.64)
ic = flopy.mf6.ModflowGwfic(gwf, strt=15.0)


Von hieran werden die Randbedingungen definiert. Eine solche ist z.B die Dirichlet-Randbedingung, die Zellen einen festen Wert zuteilt. Stell dir einen unendlichen Wasserzufluss vor, der von außen deinem Modellgebiet zufließt. Ein Pumpen in der Mitte ändert die Wasserstände am Modellrand nicht.
Die Zellen am Modellrand finden wir mit einer for-Schleife, die eine if-Abrage enthält. Sie geht die rows und columns einmal von oben nach unten und von links nach rechts durch und checkt, ob es sich um Randzellen handelt. Diese werden dann als Tupel ([(layer, row, col), head_wert]) in eine Liste gepackt (chdata) und ihnen ein fester Wasserstand zugewiesen. Schließlich erstellen wir mit der FloPy-Funktion "ModflowGwfchd" das CHD-Paket, das MODFLOW erwartet.   Überspringt man diesen Teil, behandelt MODFLOW alle Randbereiche als No-Flow-Boundary, sprich die Badewanne liefe leer irgendwann.
 Im ausklappbaren Teil unten ist eine Alternative Version, die oben einen höheren und im unteren Randbereich einen niedrigeren festen Wasserstand hat sowie links und rechts eine No-Flow-Boundary.
import os

import numpy as np
import matplotlib.pyplot as plt
import flopy

ws = "./baugrube_model"                            
sim_name = "baugrube"                             
exe_path = os.path.expanduser(r"PFAD/modflow_bin/mf6")

sim = flopy.mf6.MFSimulation(sim_ws=ws, exe_name=exe_path)
tdis = flopy.mf6.ModflowTdis(sim, nper=1, perioddata=[(1.0, 1, 1.0)])
ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE")
gwf = flopy.mf6.ModflowGwf(sim, modelname=sim_name)

  nlay, nrow, ncol = 1, 50, 50
dis = flopy.mf6.ModflowGwfdis(
    gwf,
    nlay=nlay, nrow=nrow, ncol=ncol,
    delr=10.0, delc=10.0,
    top=20.0, botm=0.0          
)

npf = flopy.mf6.ModflowGwfnpf(gwf, k=8.64)
ic = flopy.mf6.ModflowGwfic(gwf, strt=15.0)


  chd_data = []
for r in range(nrow):
  for c in range(ncol):
        if r == 0 or r == nrow - 1 or c == 0 or c == ncol - 1:
            chd_data.append([(0, r, c), 15.0])  

chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chd_data)

Alternative
h_top = 18.0   
h_bottom = 12.0 

chd_data = []

for c in range(ncol):
    chd_data.append([(0, 0, c), h_top])
    
for c in range(ncol):
    chd_data.append([(0, nrow - 1, c), h_bottom])


chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chd_data)
  



Als nächstes definieren wir eine Variable für die Pumprate und die Lokation der Brunnen in unserem Gitter. Das WEL-Paket von MODFLOW erwartet dabei folgende Datenstruktur [(layer row column),pumping_rate], die wir hier für alle Brunnen als Liste "wel_data" anlegen. Mit der Funktion "ModflowGwfwel" erstellen wir schließlich das WEL-Paket für den MODFLOW-Kernel.


import os

import numpy as np
import matplotlib.pyplot as plt
import flopy

ws = "./baugrube_model"                            
sim_name = "baugrube"                             
exe_path = os.path.expanduser(r"PFAD/modflow_bin/mf6")

sim = flopy.mf6.MFSimulation(sim_ws=ws, exe_name=exe_path)
tdis = flopy.mf6.ModflowTdis(sim, nper=1, perioddata=[(1.0, 1, 1.0)])
ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE")
gwf = flopy.mf6.ModflowGwf(sim, modelname=sim_name)

  nlay, nrow, ncol = 1, 50, 50
dis = flopy.mf6.ModflowGwfdis(
    gwf,
    nlay=nlay, nrow=nrow, ncol=ncol,
    delr=10.0, delc=10.0,
    top=20.0, botm=0.0          
)

npf = flopy.mf6.ModflowGwfnpf(gwf, k=8.64)
ic = flopy.mf6.ModflowGwfic(gwf, strt=15.0)

chd_data = []
for r in range(nrow):
    for c in range(ncol):
        if r == 0 or r == nrow - 1 or c == 0 or c == ncol - 1:
            chd_data.append([(0, r, c), 15.0])  

chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chd_data)

pumping_rate = -450.0  

wel_data = [
    [(0, 20, 20), pumping_rate],
    [(0, 20, 29), pumping_rate], 
    [(0, 29, 20), pumping_rate],  
    [(0, 29, 29), pumping_rate], 
]
wel = flopy.mf6.ModflowGwfwel(gwf, stress_period_data=wel_data)


Als letztes kommt der Output-Control-Teil (oc). Hier speichern wir die berechneten Grundwasserstände in einer .hds-Datei ab, wobei alle Grundwasserhöhen (HEAD) für alle Zeitschritte (ALL) aufgenommen werden. Mit der Funktion sim.write_simulation() übersetzen wir alle Python-Eingaben in die MODFLOW-Dateien. Mit run_simulation() führt der Code schließlich den MODFLOW-Kernel aus. Falls alles richtig läuft (konvergiert zu einer Lösung), geben wir mit dem Wahrheitswert "success" noch ne entsprechende Meldung raus.


import os

import numpy as np
import matplotlib.pyplot as plt
import flopy

ws = "./baugrube_model"                            
sim_name = "baugrube"                             
exe_path = os.path.expanduser(r"PFAD/modflow_bin/mf6")

sim = flopy.mf6.MFSimulation(sim_ws=ws, exe_name=exe_path)
tdis = flopy.mf6.ModflowTdis(sim, nper=1, perioddata=[(1.0, 1, 1.0)])
ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE")
gwf = flopy.mf6.ModflowGwf(sim, modelname=sim_name)

  nlay, nrow, ncol = 1, 50, 50
dis = flopy.mf6.ModflowGwfdis(
    gwf,
    nlay=nlay, nrow=nrow, ncol=ncol,
    delr=10.0, delc=10.0,
    top=20.0, botm=0.0          
)

npf = flopy.mf6.ModflowGwfnpf(gwf, k=8.64)
ic = flopy.mf6.ModflowGwfic(gwf, strt=15.0)

chd_data = []
for r in range(nrow):
    for c in range(ncol):
        if r == 0 or r == nrow - 1 or c == 0 or c == ncol - 1:
            chd_data.append([(0, r, c), 15.0])  

chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chd_data)

pumping_rate = -450.0  

wel_data = [
    [(0, 20, 20), pumping_rate],
    [(0, 20, 29), pumping_rate], 
    [(0, 29, 20), pumping_rate],  
    [(0, 29, 29), pumping_rate], 
]
wel = flopy.mf6.ModflowGwfwel(gwf, stress_period_data=wel_data)

oc = flopy.mf6.ModflowGwfoc(
    gwf,
    head_filerecord=f"{sim_name}.hds",
    saverecord=[("HEAD", "ALL")]
)

sim.write_simulation()
success = sim.run_simulation()

if success:
    print("Modellrechnung erfolgreich!")


Visualisieren der Ergebnisse

Wir speichern die Ergebnisse der MODFLOW-Rechnung hinter einer Variable. Erstmal definieren wir den Pfad:
dateipfad = os.path.join(ws, f"{sim_name}.hds")
ergebnis_datei = flopy.utils.HeadFile(dateipfad)

Dann laden wir alle Wasserstände aus unserem Gitter, das in unserem Fall nur einen Layer hat (deswegen 0).


dateipfad = os.path.join(ws, f"{sim_name}.hds")
ergebnis_datei = flopy.utils.HeadFile(dateipfad)

wasserstaende_3d = ergebnis_datei.get_data()
wasserstaende = wasserstaende_3d[0]

Als nächstes spannen wir eine Leinwand auf, die mit unserem Modell und den berechneten Wasserständen "bemalt" wird. Indem wir den Aspect-ratio auf "equal" setzen, verhindern wir, dass matplotlib das Gitter verzerrt. FloPy bekommt mit der Funktion "PlotMapView" die Achsen der Leinwand zur Verfügung, um das Grundwassermodell, das wir oben "gwf" genannt hatten, drauf zu malen.


dateipfad = os.path.join(ws, f"{sim_name}.hds")
ergebnis_datei = flopy.utils.HeadFile(dateipfad)

wasserstaende_3d = ergebnis_datei.get_data()
wasserstaende = wasserstaende_3d[0]

fig, ax = plt.subplots(figsize=(10, 10))
ax.set_aspect("equal") 
karten_zeichner = flopy.plot.PlotMapView(model=gwf, ax=ax)




Jetzt zeichnen wir das Gitter und die Wasserstände als Heatmap, die von Gelb über Grün zu Blau geht und die Bandbreite der Wasserstände automatisch ermittelt.
Wir lassen außerdem die Höhen berechnen und zeichnen diese als Konturlinien in Abständen 0,5m samt Beschriftung in die Abbildung.


dateipfad = os.path.join(ws, f"{sim_name}.hds")
ergebnis_datei = flopy.utils.HeadFile(dateipfad)

wasserstaende_3d = ergebnis_datei.get_data()
wasserstaende = wasserstaende_3d[0]

fig, ax = plt.subplots(figsize=(10, 10))
ax.set_aspect("equal") 
karten_zeichner = flopy.plot.PlotMapView(model=gwf, ax=ax)


karten_zeichner.plot_grid(color="lightgray", lw=0.5, alpha=0.7)

VMIN_FEST = 5.0   # z.B. minimale Ziel-Grundwasserhöhe in der Baugrube
VMAX_FEST = 15.0  # z.B. unbeeinflusster Ausgangs-Grundwasserstand

farb_flaeche = karten_zeichner.plot_array(
    wasserstaende, 
    cmap="YlGnBu_r", 
    vmin=VMIN_FEST, 
    vmax=VMAX_FEST
)

hoehen = np.arange(np.floor(wasserstaende.min()), np.ceil(wasserstaende.max()) + 0.5, 0.5)

hoehen_linien = karten_zeichner.contour_array(wasserstaende, levels=hoehen, colors="black", linewidths=0.8)
plt.clabel(hoehen_linien, inline=True, fontsize=8, fmt="%.1f m")



Schließlich zeichnen wir die Baugrube und die Brunnen und fügen eine Legende sowie Achsenbeschriftungen hinzu. Dazu definieren wir erstmal die Umrisse der Baugrube und die Lokationen der Pumpen. Zu guter Letzt lassen wir uns noch den Wasserstand in der Baugrubenmitte anzeigen.
dateipfad = os.path.join(ws, f"{sim_name}.hds")
ergebnis_datei = flopy.utils.HeadFile(dateipfad)

wasserstaende_3d = ergebnis_datei.get_data()
wasserstaende = wasserstaende_3d[0]

fig, ax = plt.subplots(figsize=(10, 10))
ax.set_aspect("equal") 
karten_zeichner = flopy.plot.PlotMapView(model=gwf, ax=ax)


karten_zeichner.plot_grid(color="lightgray", lw=0.5, alpha=0.7)

VMIN_FEST = 5.0   # z.B. minimale Ziel-Grundwasserhöhe in der Baugrube
VMAX_FEST = 15.0  # z.B. unbeeinflusster Ausgangs-Grundwasserstand

farb_flaeche = karten_zeichner.plot_array(
    wasserstaende, 
    cmap="YlGnBu_r", 
    vmin=VMIN_FEST, 
    vmax=VMAX_FEST
)

hoehen = np.arange(np.floor(wasserstaende.min()), np.ceil(wasserstaende.max()) + 0.5, 0.5)

hoehen_linien = karten_zeichner.contour_array(wasserstaende, levels=hoehen, colors="black", linewidths=0.8)
plt.clabel(hoehen_linien, inline=True, fontsize=8, fmt="%.1f m")


baugrube_x = [200, 300, 300, 200, 200]
baugrube_y = [200, 200, 300, 300, 200]
ax.plot(baugrube_x, baugrube_y, color="red", linestyle="--", linewidth=2, label="Grundwasserabsenkung")

brunnen_x = [205, 295, 205, 295]
brunnen_y = [205, 205, 295, 295]
ax.plot(brunnen_x, brunnen_y, "r^", markersize=10, label="Absenkbrunnen")


plt.colorbar(farb_flaeche, shrink=0.8, label="Grundwasserstand [m]")

plt.title("Baugrubenwasserhaltung")
plt.xlabel("X-Koordinate [m]")
plt.ylabel("Y-Koordinate [m]")
plt.legend(loc="upper right")

plt.show() 

wasserstand_mitte = wasserstaende[25, 25]
print(f"Wasserstand in der Baugrubenmitte: {wasserstand_mitte:.2f} m")

Wir brauchen einen Zusammenhang der erklärt, wovon die nötige Förderrate abhängig ist. Also eine Funktion, die so aussehen muss: Q= Irgendwelche Parameter, die sich gegenseitig beeinflussen. 







import numpy as np # brauchen wir um Pi und eine range von Zahlen zu generieren. 

import matplotlib.pyplot as plt # brauchen wir um Diagramme zu plotten.

radien = np.arange(1, 11)      # erstellt Zahlenreihe von 1 bis 10, in diesem Fall Zylinderradien
hoehen = np.arange(1, 11)      # # erstellt Zahlenreihe von 1 bis 10, in diesem Fall Zylinderhöhen

mantelflaeche = 2 * np.pi * radien * hoehen

Plot
plt.plot(radien, mantelflaeche, marker='o')
plt.title("Mantelfläche in Abhängigkeit des Radius und der Höhe")

plt.xlabel("Radius [m]")
plt.ylabel("Mantelfläche [m²]")
plt.grid(True)
plt.show



Das heißt je größer die Fläche ist, die einen gewünschten Absenkungsziel entsprechen soll, desto mehr muss abgepumpt werden. Das kann man ausdrücken durch Q=(2*pi*r*h*K*(H-h))/ln R. 

Grundwasserdruckhöhe im gespannten Grundwasserleiter.

Bei einem vollverfilterten Brunnen wird ein mit der Zeit immer kleinerer Teil der Filterstrecke durchströmt, sodass der
  
Durchflussquerschnitt ändert sich bei Grundwasser mit gespannter Oberfläche nicht mit der Entfernung zum Brunnen. 
  
 
  
  Brunnenradius. was ist der unwirksame Brunnenradius?
  




Warum natürlicher Logarithmus der Reichweite und des wirksamen Brunnenradius?
Was drückt der Logarithmus aus? Wann braucht man ihn mathematisch?

. Tabelle plotten. Mathematisch durch lograithmische Funktion beschrieben.

Dann bleibt die Frage, wie weit wir absenken wollen? 


Überlegung. ganz viele Zylinder um den Brunnen zeichnen und schauen wie sich was ändert.

logarithmic functions. These come in handy when we need to consider any phenomenon that varies over a wide range of values, such as pH in
chemistry or decibels in sound levels oder eben die Reichweite eines Absenkungstrichter. Eine Funktion, die mit der Entfernung  proportional zur aktuellen Entfernung abnimmt ist eine Exponentialfunktion. Heißt?