Excel hat Zellen, denen man einen Anfangswert und eine Beziehung zu den benachbarten Zellen zuweisen kann. 
        
                             Grundwasserneubildung
        Inflow Linke Zelle ->Aktive Zelle ->Rechte Zelle Outflow
                             Brunnen
        Change in Groundwater storage ist, was in der mittleren Zelle übrig bleibt. Was ist das?
        
      

      Virtuelle Umgebung erstellen: Man erstellt einen Ordner, wo man alles (Python-Kernel, Bibliotheken, Module) reinpackt. Dazu steuern wir im Terminal mit "cd" den Pfad an und erstellen dort mit einem Modul (-m) von Python 3.13 eine virtuelle Umgebung genannten Ordner, hier genannt "virtuelleumgebung". Durch das Aktivieren wird der Python-Kernel nur in diesem Ordner angesteuert. Mit pip installieren wir schonmal alle Bibliotheken, die später nützlich sein werden. Das geht auch später noch im jupyter notebook, indem man ein "!pip install" schreibt, also mit Ausrufezeichen davor. registrieren den Kernel für jupyter notebooks
    
       
orkan@linux:~> cd /home/orkan/Downloads/
orkan@linux:~/Downloads> python3.13 -m venv virtuelleumgebung
orkan@linux:~/Downloads> source virtuelleumgebung/bin/activate
(virtuelleumgebung) orkan@linux:~/Downloads> pip install scipy numpy matplotlib ipykernel notebook flopy
(virtuelleumgebung) orkan@linux:~/Downloads> python -m ipykernel install --user --name="flopy_env" --display-name "Python (FloPy)"
(virtuelleumgebung) orkan@linux:~/Downloads> jupyter notebook

    
    Download latest MODFLOW und exportiere die Datei in einen Pfad deiner Wahl.
    
    Why I prefer Open Source or Commons? To move money from licences to Human Capital thats why.
    

      
      
      
      
      "File" - "Import" - "Shapefile..."
      "Object" - "Show or Hide Objects..." - Right-click on Shape-File - "Edit" - Right-click on "Unused Objects" - "Object Properties" - "Color object line" - Toggle on "Use to Set Grid Size" - "Data Sets" - Define your Shapefile as as Model_Top.
      
      "Data" - "Edit Data Sets" - "Required" - "Hydrology" - "Modflow_Initial_Head"
      
      Like Marco Pierre White said: "I do it this way. I don't know why. It actually makes it harder." 
      
      
      

1
      10 km breites Tal und viele Kilometer lang. 3 Schichten. Oben ungespannter Grundwasserleiter aus Sand. Dann eine schluffiger Aquitard, der den unteren sandigen gespannten Grundwasserleiter von dem obigen trennt. Darunter ist undurchlässiges Festgestein. Im Westen und Osten undurchlässige Hügel. Im Osten fließt ein Fluss, der an den obigen Grundwasserleiter über seine ganze Mächtigkeit angrenzt. Er habe über die gesamte Länge einen Pegelstand von 320m, also eine Constant-Head-Boundary. Ein Kanal mit einem Pegelstand von 330m begrenzt auf gleiche Weise den oberen Grundwasserleiter im Westen. Es wird angemerkt, dass diese Annahme das Modell vereinfacht (ich glaube weil dann nur horizontaler Strom herrscht und nicht noch von unten der Fluss durchströmt wird.) und "Such deep waterways are not typical of the real world". Im Norden und Süden liegen kilometerentfernt Moränen, die undurchlässig sind. Weil sie orthogonal zur Strömungsrichtung liegen, sind es No-Flow-Boundaries. Hydraulische Leitfähigkeiten werden angegeben: Oben in horizontaler Richtung 50 m/d, in vertikaler Richtung 10 m/d, Speicherkoeffizient 0,2. Basis des oberen Aquifers liegt auf 220m, Top entspricht der Geländeoberfläche bei 350m. Aquitard hat ne hydraulische Leitfähigkeit von 0,01 m/d. Speicherkoeffizient 10^-6. Basis des Aquitards liegt bei 200 m. Unterer Grundwasserleiter hat ne hydraulische Durchlässigkeit von 200 m/d in horizontaler Richtung und 20 m/d in vertikaler Richtung. Speicherkoeffizient liegt bei 10^-6. Die Basis dieses Grundwasserleiters sei bei 0m.
      
      Aufgabe ist, in einer Area-of-Interest (eine quadratische Farm mit Seitenlänge 4 km genau zwischen Fluss und Kanal gelegen) die Grundwasserströmung im Gleichgewichtszustand zu modellieren.
      
      Use 3 layers, 21 rows and 20 columns with a uniform horizontal grid spacing of 500 m each direction.
      
      
      
      
      Wie immer importieren wir zuerst alle Module, die wir brauchen werden:
      import os
import matplotlib.pyplot as plt
import numpy as np
import flopy

Wir legen nun fest, wo die Input- und Output-Dateien für MODFLOW (z. B. .nam, .dis, .hds) gespeichert werden sollen, welchen Namen sie tragen sollen und wo der MODFLOW-Rechenkern steckt aufm Computer.
ws = "./übung3"
sim_name = "übung3"
exe_path = os.path.expanduser("~/modflow_bin/mf6")

Dann erstellen wir per flopy-Befehle, die MODFLOW-Simulationsdatei, die alle Pakete miteinander verknüpft, die MODFLOW braucht.
Das Paket für Zeitdiskretisierung (TDIS) enthält Infos über die Anzahl der Stressperioden (hier nper=1), Länge der Stressperiode (perlen=1.0), Unterteilung der Stressperiode (nstp = 1) und dann noch ein time-step Multiplikator für Zeitschritte (hier sind mit tsmult= 1.0 alle Zeitschritte gleich lang).
Schließlich erstellen wir noch den Gleichungslöser IMS (Iterative Model Solver) und wählen hier SIMPLE für bestimmte Konvergenzkriterien und Iterationszahlen.
Letztlich noch das Groundwater-Flow-Paket, an das wir die Geometrie später ranhängen. Ja, es ist alles verdammt abstrakt, aber man muss da durch.
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)

Jetzt bauen wir die Geometrie auf. Wir definieren erstmal Anzahl der Layer, Reihen und Spalten. Dann eine Liste mit den Unterkanten der drei Layer. Erst dann erstellen wir das Diskretisierungspaket und hängen es an das Groundwater-Flow-Paket. Die erzeugten Quadrate sollen 500 m Seitenlänge haben.
nlay, nrow, ncol = 3, 21, 20
botm_layers = [220, 200, 0.0]
dis = flopy.mf6.ModflowGwfdis(
    gwf,
    nlay=nlay,
    nrow=nrow,
    ncol=ncol,
    delr=500,
    delc=500,
    top=350.0,  
    botm=0.0,
)

Jetzt definieren wir die hydraulischen Durchlässigkeiten für jeden Layer in horizontaler und vertikaler Richtung und deklarieren, ob gespannt (1; durchströmte Mächtigkeit bleibt konstant) oder ungespannt (0; Transmissivität ändert sich). Dann erstellen wir das Node-Property-Flow-Paket, das die eben definierten Eigenschaften dem Modell zuweist.

k_horizontal = [50.0, 0.01, 200.0]
k_vertical = [10.0, 0.01, 20.0]

icelltype_layers = [1, 0, 0]

npf = flopy.mf6.ModflowGwfnpf(
    gwf, k=k_horizontal, k33=k_vertical, icelltype=icelltype_layers
)

Mit dem Initial-Conditions-Paket geben wir einen Anfangswasserstand vor (Hier einen Wert, der zwischen den Oberflächengewsserpegeln liegt). Und jetzt kommt der eklige Teil, wenn man ohne GUI arbeitet: Randbedingungen definieren in den entsprechenden Zellen. Dadurch, dass es sich um ne symmetrische Geschichte handelt, ist das hier noch ok. Wir machen ne Liste mit den Zellen, die einen fixen Wasserstand bekommen sollen. Dazu gehen wir alle Zellreihen durch, die im oberen Layer in der Spalte 0 (westliche Grenze) und in der letzten Spalte sind (östliche Grenze) und weisen ihnen entsprechende Wasserstände zu. Damit erstellen wir das CHD-Paket. Die letzte Spalte hat den Index 19, weil python bei 0 beginnt zu zählen. Weil wir oben ja ncol=20 definiert hatten ergibt sich hier ncol-1. Die Nord-, Süd- und Sohlränder sind nicht definiert und damit automatisch No-Flow-Boundaries.

ic = flopy.mf6.ModflowGwfic(gwf, strt=325.0)

chd_data = []

for r in range(nrow):
    
    chd_data.append([(0, r, 0), 330.0])

    # Fluss am östlichen Rand (Spalte 19, Schicht 0): Wasserstand 320 m
    chd_data.append([(0, r, ncol - 1), 320.0])

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

Als letztes erstellen wir ein Paket(Output Control Package), das die Ausgabedaten koordiniert. Was auch immer das heißt.
Wir geben eine Binärdateifür die berechneten Grundwasserstände aus, eine Binärdatei für die Wasserbilanz und Volumenflüsse von Zelle zu Zelle (cell by cell-Flux). Außerdem speichern wir alle berechneten Potentiale und Bilanzdaten. Schließlich übersetzen wir unsere Eingabedaten in ASCII-Dateien, die MODFLOW versteht und starten die Berechnung, die, wenn sie glückt, ein True ausgibt.

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

sim.write_simulation()
success = sim.run_simulation()
Die MODFLOW listing-file zeigt die Ergebnisse der Berechnung in einer .lst-Datei, die als .txt-Datei gelesen werden kann. Sie findet sich im weiter oben definierten Ablageordner. Wir können dort unter "VOLUME BUDGET FOR ENTIRE MODEL AT END OF TIME STEP" sehen, wie viel Wasser aus dem Kanal strömt, undzwar: 96990 m^3/Tag über die gesamte Strecke von 21 Reihen x 500 m Zellgröße, also 10500 m. Wir sehen auch, das genau diese Menge in den Fluss geht, wir also ne "PERCENT DISCREPANCY" von 0 haben.





Berechnung des Einzugsgebiets eines Brunnens

Übung aus ISBN 978-3-443-01032-4

West-Ost-Erstreckung des Aquifers 2500 m. Westgrenze No-Flow-Boundary. Ostgrenze ein Fluss in vollständigem hydraulischem Kontakt mit dem Aquifer und einem festen Wasserstand, also ein Festpotentialrand. Nach Süden und Norden hin sind es mehrere Kilometer. Aquifer wird in x-Richtung mit 50 und in y-Richtung mit 51 Zellen angegeben. Zellabmessungen in jede Richtung 50 m. Aquifer näherungsweise gespannt. Homogen und isotrop. Transmissivität = 0.005 m2/s. Aquifermächtigkeit = 10 m.  Grundwasserneubildung wird mit 8 L/s*km^2 (=8*10^-9 m^3/s*m^2) angegeben. Entnahmerate eines Brunnens -0.005 m^3/s. Brunnen ist vom Westrand 1450 m und vom Vorfluter 1000 m entfernt. Vorfluter hat auf Höhe des Brunnens eine Breite von 50 m. Durchflusswirksame Porosität nf=0.1 gilt für das gesamte Gebiet. Frage. Wie groß ist das Einzugsgebiet des Brunnens bei stationären Verhältnissen? Einmal numerisch und dann analytisch.

  Bauen wir es einmal in flopy nach. Wie jedes gute Gericht angefangen mit Knoblauch, Zwiebeln und Olivenöl:
import os
import matplotlib.pyplot as plt
import numpy as np
import flopy

Hier taufen wir unser Baby:
ws = "./übung6.17"
sim_name = "übung6.17"
exe_path = os.path.expanduser(r"C:\WRDAPP\bin\mf6")

Das Modell soll stationär sein, also einen zeitunabhängigen Gleichgewichtszustand darstellen bei dem sich die Randbedingungen (z. B. Grundwasserneubildung, Flusswasserstände) über die Zeit nicht ändern. Das heißt, wir nehmen eine Stressperiode (nper=1) und nur einen Zeitschritt (nstp = 1).
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)

Wir übertragen die Geometrie des Grundwasserleiters:

nrow, ncol = 50, 50
delr, delc = 50.0, 50.0
top, botm = 10.0, 0.0

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

In der Aufgabe ist die Transmissivität angegeben. Die Transmissivität ist das Produkt aus hydraulischer Durchlässigkeit und Aquifermächtigkeit. Weil wir gespannte Verhältnisse haben, ist letztere konstant, sodass wir auf den k-Wert schließen können.

npf = flopy.mf6.ModflowGwfnpf(
    gwf,
    k=0.005 / 10.0,
    icelltype=0
)

Wir legen die Randbedingungen fest im Format: [(layer, row, col), head_value]. Jede Grenze, der keine explizite Randbedingung (wie CHD, RIV, WEL, GHB) zugewiesen wird, ist automatisch ein No-Flow Boundary. Zuerst den Anfangswasserstand.

ic = flopy.mf6.ModflowGwfic(gwf, strt=10.0)

# Ostgrenze: 50 m breiter Fluss mit festem Wasserstand (CHD = 10.0 m)
chd_data = [[(0, r, ncol - 1), 10.0] for r in range(nrow)]

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


Damit Grundwasser strömen kann, brauchen wir irgendwo eine Wasserstandsdifferenz, also einen Wasserein- und/oder Austrag. In unserem Beispiel kommt Wasser durch Regen und geht aus einer Pumpe:

  rch = flopy.mf6.ModflowGwfrcha(gwf, recharge=8e-9)

wel_data = [[(0, 25, 29), -0.005]]

wel = flopy.mf6.ModflowGwfwel(gwf, stress_period_data=wel_data)

Und dann noch das übliche Speichern, ausführen

head_file = f"{sim_name}.hds"
budget_file = f"{sim_name}.cbc"

oc = flopy.mf6.ModflowGwfoc(
    gwf,
    head_filerecord=head_file,
    budget_filerecord=budget_file,
    saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
)

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

Wir wollen das ganze plotten. Dazu lesen wir zunächst die Wasserstände aus der gwf-Datei aus und speichern sie in einer Variable ("heads"). Daraus geben wir uns schon Min und Max aus.


heads = gwf.output.head().get_data()
print("Max. Wasserstand (Westen):", round(heads[0].max(), 2), "m")
print("Min. Wasserstand (Brunnen):", round(heads[0].min(), 2), "m")

 Als nächstes plotten wir uns
 
pmv = flopy.plot.PlotMapView(model=gwf)

# Konturlinien (Höhenlinien) zeichnen und mit "m" beschriften
konturen = pmv.contour_array(heads[0], colors="black")
plt.clabel(konturen, fmt="%.2f m")

Zur Orientierung zeichnen wir noch unseren Brunnen und den Fluss ein. Wir heben hervor alle Zellen, die in dem WEL- bzw. in dem CHD-Paket enthalten sind:

# Randbedingungen einzeichnen
pmv.plot_bc("WEL", color="red") 
pmv.plot_bc("CHD", color="blue") 
plt.show()


      

• Contaminant Plume Modeling with MODFLOW
• Plume Response Modeling
• Delineating Well Head Capture Zones
• Time-of-travel(TOT)-paths, Isochronen