IN-KJM1900 - H?st 2018

Koblede diffligninger

IN-KJM1900, femte forelesning, 13.november

Mål for økta:

  • Litt repetisjon om interpolasjon
  • Fortsettelse p? deloppgave 2
    • Hydrologisk delmodell
    • Kjemisk delmodell
  • Ny teori:
    • Koblede differensialligninger

Orakeltjeneste, uke 46

Orakelhjelp tilbys ved Kjemisk Institutt p? f?lgende tidspunkter denne uka

Tirsdag 6/11 14.00-16.00
Onsdag 7/11 13.00-15.00
Torsdag 8/11 13.00-15.00

Innlevering av deloppgave 2

  • Fristen for ? levere deloppgave 2 av prosjektet er s?ndag 18. november klokka 23.59.
  • Oppgaven leveres (som normalt) p? Devilry.?
    • ?
  • Du oppfordres sterkt til ? benytte gruppetimer og samretting p? torsdag klokka 10.15-12.00
  • Sp?rsm?l om prosjektet?

Interpolasjon (repetisjon)

Interpolasjon er en teknikk hvor vi benytter diskrete datapunkter for ? gj?re et mer eller mindre kvalifisert gjett p? hva som finnes mellom.

Dette er p? en m?te en "anti-diskretisering".

Interpolasjon

I dette kurset vil vi benytte en interpolasjons-funksjon fra scipy, som importeres og benyttes p? f?lgende m?te:

In [2]:
from scipy.interpolate import interp1d 
import numpy as np

# data punkter for testing
data_x = np.array([1,2,3,4])
data_y = np.array([1,7,3,5])

y_interp = interp1d(data_x, data_y, 3) #gj?r en kubisk (3. ordens) interpolasjon

x = np.linspace(1,4,100) #?k oppl?sningen p? grid
y = y_interp(x) #beregn funksjonsverdiene til den interpolerte funksjonen

Kopier gjerne koden over og kj?r p? din egen maskin for ? pr?ve den ut.

Eksamensrelevante sp?rsm?l

Prosjektets del 2: Birkenesmodellen

  • Numerisk modell for vannbevegelse og kjemi i jordv?ske.
  • Tilpasset observasjoner i nedb?rsfelt i Birkenes kommune, Aust-Agder.
  • Flere publikasjoner p? 70-, 80- og 90-tallet.
  • Kjemisk modell koblet til hydrologisk modell.

Hydrologisk delmodell: detaljer

  • Nedb?rsfelt modellert som reservoarer $\mathbf{A}$ og $\mathbf{B}$
      • Her benyttes fet font for "merkelappene" p? reservoarene.
  • For et reservoar, for eksempel merket $\mathbf{M}$:
    • Vannstand som funksjon av tid: $M(t)$
    • Str?m ut: $Q_M(t)$
    • Evaporasjon: $E_M(t)$
    • Nedb?r ("precipitation") $P_M(t)$
    • Generelt: subindeks indikerer reservoar.
  • Avrenning til bekk. Fysiske ma?linger gjort her.

Det kan l?nne seg ? skrive en liste over st?rrelsene som inng?r i problemet.

Reservoar $\mathbf{A}$

Vi beskriver f?rst modellen for reservoar A:

Vannstand som funksjon av tid: $A(t)$ m?les i $mm$, hvor $mm$ tilsvarer vannmengde per flateenhet slik at $1mm = 1 \frac{l}{m^2}$.

Vann str?mmer inn i reservoaret kun gjennom nedb?ren $P(t)$.

Vann str?mmer ut fra reservoaret gjennom evaporasjon $E_A$ og en bekk $Q_A$.

Vannstr?m/fluks m?les i $mm/dag = mm \cdot dag^{-1}$.

Reservoar $\mathbf{A}$

Regler for reservoar $\mathbf{A}$:
  • Ingen avrenning ($Q_A = 0$) n?r vannstanden $A(t)$ er mindre enn $A_{min}$
  • Avrenning intensiveres n?r vannstanden $A(t)$ ?ker over $A_{min}$
  • Evaporasjon er en funksjon av temperatur $E_A := E_A(T)$

Vi skal uttrykke dette matematisk for ? utforme modellen.

Reservoar $\mathbf{A}$

Vannstanden i $\mathbf{A}$ p?virkes positivt (?kes) av nedb?ren $P$, og negativt av evapotranspirasjon $E$ og utstr?mming til bekken $Q_A$, alts? er

$A := A(t, P_A, E_A, Q_A)$.

Vi betrakter n? en endring i vannstanden over (i l?pet av) et lite tidsintervall $\delta t$ i magasin $\mathbf{A}$, angitt ved $\delta A$:

$\frac{\delta A}{\delta t} = P - E_A - Q_A \leftrightarrow \delta A = P \delta t - E_A \delta t - Q_A \delta t$

Legg merke til:

  • Symbolet $\delta$ er en liten gresk delta ($\Delta$), og brukes konvensjonelt for ? beskrive endringer og differanser. Andre steder kan vi bruke $\Delta$ eller noen ganger $h$.
  • Dimensjonen til $\delta t$ er tid, vannstr?m/fluks har dimensjon vannmengde per flateenhet per tid ($mm/d?gn$) og $\delta A$ har dimensjonen vannmengde per flateenhet.
  • Sp?rsm?l til diskusjon: g?r dimensjonene i uttrykket for endringen opp? Diskuter.

$$ mm = \frac{mm}{d?gn} \cdot \text{d?gn} + \frac{mm}{d?gn} \cdot \text{d?gn} - \frac{mm}{d?gn} \cdot \text{d?gn} $$

Str?m $Q_A$ som funksjon av vannstanden $A(t)$

N?r vannstanden $A(t)$ er h?y er str?mmen ut av $\mathbf{A}$ sterk, og n?r vannstanden er lav er str?mmen svak.

Str?mmen er derfor en funksjon av vannstanden $A(t)$, og i modellen v?r uttrykkes dette som

$Q_A := K_A \text{max}(A - A_{min}, 0)$,

hvor $K_A = 0.8 \text{d?gn}^{-1}$ og $A_{min} = 13mm$.

Evapotranspirasjon $E_A$ som funksjon av temperaturen $T$

Ved h?ye temperaturer g?r mange prosesser fortere.

I dette tilfellet har vi h?yere fordamping og ?kt opptak i vegetasjon, og vi uttrykker det som en line?r relasjon:

\begin{equation} E_{A}(A, T) = \left\{?\begin{array} 00.2\cdot T &~~& A > 1mm \\ 0 & & \mbox{ellers} \end{array} \right. \end{equation}

Legg merke til at vi ogs? har en avhengighet av $A(t)$ her.

En differensialligning for reservoar $\mathbf{A}$

Igjen, med utgangspunkt i

$\frac{\delta A}{\delta t} = P(t) - E_A(A,t) - Q_A(A,t) $

St?rrelsen $\frac{\delta A}{\delta t}$ kan leses som "en liten endring i $A$ som f?lge av en liten endring i $t$". Om vi g?r fra "en liten endring" til en "infinitesimal endring" f?r vi en differensialligning:

$\frac{d}{dt}A(t) = P(t) - E_A(A, t) - Q_A(A, t) $

En differensialligning for reservoar $\mathbf{A}$

$\frac{d}{dt}A(t) = P(t) - E_A(A, t) - Q_A(A, t) $

Kan vi l?se ligningen analytisk?

Nei, blant annet fordi nedb?ren

$P(t)$

(og temperatur) skal leses inn fra en tabell. Ligningen m? derfor l?ses numerisk som et initialverdiproblem.

En numerisk l?sning av diffligningen for reservoar $\mathbf{A}$

Vi begynner som vanlig med ? diskretisere variablene, i dette tilfellet tiden

$$t \rightarrow t_i = i \cdot \Delta t$$

Vi benytter deretter en numerisk approksimasjon p? den tidsderiverte (fremoverdifferansen):

$$\frac{d}{dt}f(t) \approx D^+f(t) = \frac{f(t + \Delta t) - f(t)}{\Delta t}$$

Dette gir

$\frac{d}{dt}A(t_n) \approx D^+A(t_n) = \frac{A(t_{n+1})-A(t_n)}{\Delta t} = P(t_n) - E_A(t_n) - Q_A(t_n) $

Tenk over: Hva er kjent/ukjent i denne ligninga? Hva er det vi vil l?se for?

En numerisk l?sning av diffligningen for reservoar $\mathbf{A}$

Til slutt l?ser vi den ukjente $A(t_{n+1})$ algebraisk og finner:

$$A(t_{n+1}) = A(t_n) + P(t_n)\Delta t - E_A(t_n)\Delta t - Q_A(t_n)\Delta t $$

Eller for ? ikke skrive for mye

$$A_{n+1} = A_n + P_n\Delta t - E_{A,n}\Delta t - Q_{A,n}\Delta t $$

Dette kalles et oppdateringsskjema, som vi ogs? husker fra forrige gang.

$\rightarrow$ Live kodeeksempel

Reservoar $\mathbf{B}$

Vi beskriver n? modellen for reservoar $\mathbf{B}$:

Vannstand som funksjon av tid: $B(t)$ m?les i $mm$.

Vann str?mmer inn i reservoaret kun gjennom bekken $Q_A(t)$.

Vann str?mmer ut fra reservoaret gjennom evaporasjon $E_B$ og en bekk $Q_B$.

Reservoar $\mathbf{B}$

Regler for reservoar $\mathbf{B}$:
  • Innstr?mming fra $Q_A$ avtar n?r vannstanden i $\mathbf{B}$ tiltar. ($A_{sig}$ er en funksjon av $B(t)$ )
  • Ingen utstr?mming til bekk n?r vannstanden er under $B_{min}$.
  • $\mathbf{B}$ har i tillegg en maksgrense $B_{max}$, hvor alt overskytende vann g?r direkte i bekken $Q_{over}$
  • Evapotranspirasjon er en funksjon av temperatur $E_B := E_B(T)$

Reservoar $\mathbf{B}$

Vannstanden i $\mathbf{B}$ p?virkes positivt (?kes) av innstr?mming fra $\mathbf{A}$ (gjennom bekken $Q_A$, og negativt av evaporasjon $E_B$ og utstr?mming til bekken $Q_B$ og $Q_{over}$. Alts? er

$B := B(t, P_A, E_A, Q_A, Q_{over})$.

Vi kan uttrykke en endring i vannstanden i $\mathbf{B}$ i l?pet av et lite tidsintervall $\delta t$

$\frac{\delta B}{\delta t} = A_{sig}Q_A - E_B - Q_B - Q_{over}$

Legg merke til fortegnene.

To viktige forskjeller fra uttrykket for $A(t)$:

  • $A_{sig}$
  • $Q_{over}$

Parameteren $A_{sig}$ er en funksjon av $B(t)$, $B_{min}$ og $B_{max}$:

\begin{equation} A_{sig} = \left\{\begin{array} 1 & & B \leq B_\text{min} \\ 1- 0.25\cdot\frac{B-B_\text{min}}{B_\text{max}-B_\text{min}} & & B_\text{min} < B \leq B_\text{max} \\ 0 & & B_\text{max} < B \end{array} \right. \end{equation}

Utstr?mming fra $\mathbf{B}$: $Q_B$

Som for $\mathbf{A}$ er ogs? str?mmen $Q_B$ ut fra $\mathbf{B}$ en funksjon av vannstanden:

$Q_B := K_B \text{max}(B - B_{min}, 0)$,

hvor $K_B = 0.045 \text{d?gn}^{-1}$ og $B_{min} = 40mm$.

Utstr?mming fra $\mathbf{B}$: $Q_{over}$

I tillegg vil, dersom vannstanden overskrider et maksniv? $B_{max}$, alt overskytende vann g? rett over i bekken:

$Q_{over} = \text{max}\left( (B - B_{max})\delta t^{-1} + A_{sig}Q_A - Q_B - E_B, 0 \right)$,

Du kan enten f?rst beregne $B(t)$ og deretter la overskytende vann g? tilbake i bekken, eller du kan benytte alternative uttrykk oppgitt i oppgaven.

Differensialligning for $B(t)$

Fra uttrykket
$\frac{\delta B}{\delta t} = A_{sig}Q_A - E_B - Q_B - Q_{over}$
finner vi igjen en differensialligning

$\frac{d}{dt} B(t) = A_{sig}(B) Q_A(A) - E_B - Q_B(B) - Q_{over}(B) $

Hva skjer her?

Vi ser p? de to differensialligningene vi har funnet for modellen s? langt:

$\frac{d}{dt}A(t) = P(t) - E_A(A) - Q_A(A) $

$\frac{d}{dt} B(t) = A_{sig}(B) Q_A(A) - E_B - Q_B(B) - Q_{over}(B) $

Koblede ligninger

Vi ser p? de to differensialligningene vi har funnet for modellen s? langt:

$\frac{d}{dt}A(t) = P(t) - E_A - \color{blue}{Q_A(A)} $

$\frac{d}{dt} B(t) = A_{sig}(B) \color{blue}{Q_A(A)} - E_B - Q_B(B) - Q_{over}(B) $

Legg merke til at vi ikke kan l?se for $B$ uten ? kjenne til $A$. Dette kalles koblede differensialligninger, eller systemer av ligninger.

Numerisk l?sning av systemet

Vi kan l?se numerisk for $B(t)$ p? samme m?te som for $A(t)$, og finner at ligningssettet v?rt er

\begin{align} A_{n+1} &= A_n + \cdot P_n \Delta t - \cdot E_{A,n}\Delta t - \cdot Q_{A,n} \Delta t. \\ B_{n+1} &= B_n + \cdot A_{\text{sig},n} \cdot Q_{A,n} \Delta t - \cdot E_{B,n}\Delta t - \cdot Q_{B,n} \Delta t - \cdot Q_{\text{over},n} \Delta t . \end{align}

Numerisk l?sning av ligningssettet

Hvilke st?rrelser kjenner vi ved f?rste iterasjon?

\begin{align} A_{1} &= A_0 + P_0 \Delta t - E_{A,0}\Delta t - Q_{A,0} \Delta t. \\ B_{1} &= B_0 + A_{\text{sig},0} Q_{A,0} \Delta t - E_{B,0}\Delta t - Q_{B,0} \Delta t - Q_{\text{over},0} \Delta t . \end{align}

Str?m i bekken

Ved ? l?se systemet iterativt (med Eulers metode) bestemmes blant annet $A(t)$, $B(t)$, $Q_A(t)$, $Q_B(t)$ og $Q_{over}(t)$.

Str?mmene inn og ut av bekken kan benyttes til ? bestemme den totale str?mmen i bekken:

$$ Q(t) = (1-A_\text{sig}(t))\cdot Q_{A}(t) + Q_{B}(t) + Q_\text{over}(t)$$.

St?rrelsen $Q(t)$ kan sammenlignes med data i kolonnen "avrenn" i birkenes.dat-filen.

Kjemisk delmodell

Viktige antagelser i den kjemiske delmodellen:

  • Sulfat kommer med nedb?ren og f?lger vannstr?mmen uten ? reagere med omgivelsene.
  • Unntak: sulfat forlater ikke systemet gjennom evapotranspirasjon.
  • Elektron?ytralitet: summen av ladningene i vannet er lik 0.
  • Likevekt: for andre ioner enn sulfat antar vi at de er fullstendig bestemt av likevektsreaksjoner.

Sulfattransport

I kolonne 1 i birkenes.dat oppgis sulfatkonsentrasjonen i nedb?ren:

  • Konsentrasjonen varierer fra dag til dag
  • Av tekniske ?rsaker er den satt til negativ verdi n?r nedb?ren P = 0.
    • Tenk over: kan dette bli et problem?
  • Konsentrasjon $C$ angis med enheten $mol \cdot l^{-1}$ (mol per liter).
  • Mengden $M$ angis (mot normalt) som $mol \cdot m^{-2}$ (mol per kvadratmeter).

Differensialligninger for sulfattransport

Ettersom evapotranspirasjon ikke transporterer sulfat ut av systemet, finner vi differensialligningene

\begin{align} \frac{d}{dt} M_A & = P\cdot C_P - Q_{A}\cdot C_A, \\ \frac{d}{dt} M_B & = A_\text{sig}\cdot Q_{A}\cdot C_A - (Q_{B} + Q_\text{over})\cdot C_B, \end{align}

St?rrelsene $Q_A$, $Q_B$, $Q_{over}$ og $A_{sig}$ har vi l?st for i den hydrologiske modellen, s? de resterende st?rrelsene kan l?ses numerisk.

G?r enhetene opp?

Vi ser p? enhetene i ligningen:

\begin{align} \frac{d}{dt} M_A & = P\cdot C_P - Q_{A}\cdot C_A \end{align}

Ja:

\begin{align} \rightarrow \frac{mol}{d?gn} = \frac{l}{d?gn} \frac{mol}{l} - \frac{l}{d?gn} \frac{mol}{l} \end{align}

Den andre ligningen f?lger samme m?nster.

Differensialligninger for sulfattransport

Vi kan l?se ligningene for sulfattransport med Eulers metode, og finner da konsentrasjon og stoffmengde sulfat i reservoarene $\mathbf{A}$ og $\mathbf{B}$ og dessuten i avrenningen $Q$.

Utfordringen er n? ? bestemme konsentrasjonen av den andre ionene.

Elektron?ytralitet

I alle vannmengder i modellen gjelder prinsippet om elektron?ytralitet:

\begin{equation} \label{elnoy} {[H^{+}] + 2[Ca^{2+}] + 3[Al^{3+}] = 2[SO^{2-}_{4}] + [HCO^{-}_{3}]} \end{equation}

Her kjenner vi ved f?rste ?yekast kun sulfatkonsentrasjonen $[SO^{2-}_{4}]$.

Hvordan finnes de andre?

Likevektrelasjoner:

\begin{align} K_{HCa} & = \frac{[{H^+}]^{2}}{[{Ca^{2+}}]} \\ K_{AlH} & = \frac{[{Al^{3+}}]}{[{H^+}]^{3}} \label{eq2} \\ K_{H} & = [{H^+}][{HCO^{-}_{3}}]\label{eq3} \end{align}

Elektron?ytralitet

Fra likevektene kan vi l?se for $[H^+]$ og finne

\begin{align} [{H^+}] + 2\frac{[{H^+}]^2}{K_{HCa}} + 3K_{AlH}[{H^+}]^3 = 2[{SO^{2-}_{4}}] + \frac{K_{H}}{[{H^+}]}\label{elnoy2} \end{align}

Vi lar $[{H^+}] := x$ og $[{SO^{2-}_{4}}] := s$ slik at uttrykket kan skrives

$$3K_{AlH}x^3 + 2K_{HCa}^{-1}x^2 + x - 2s - K_{H} x^{-1} = 0$$

Elektron?ytralitet

For et gitt tidssteg og reservoar kjenner vi alle st?rrelsene med unntak av $x$:

\begin{align} 3K_{AlH}x^3 + 2K_{HCa}^{-1}x^2 + x - 2s - K_{H} x^{-1} = 0 \end{align}

Ligningen kan derfor l?ses med Newtons metode.

Ioner i vannmassene

Om vi kan bestemme $[{H^+}]$, s? kan vi benytte likevektene til ? finne konsentrasjonen av de andre ionene:

\begin{align} K_{HCa} & = \frac{[{H^+}]^{2}}{[{Ca^{2+}}]} \\ K_{AlH} & = \frac{[{Al^{3+}}]}{[{H^+}]^{3}} \\ K_{H} & = [{H^+}][{HCO^{-}_{3}}] \end{align}

Litt om neste uke

Vi har n? gjennomg?tt de mest sentrale delene av pensum.

Det gjenst?r n? ? sette det inn i en helhetlig sammenheng, blant annet ved ? betrakte

  • Newtons metode og Eulers metode som rekurrensrelasjoner
  • Eulers metode som numerisk integrasjon

Vi vil bruke den siste forelesningen til dette, repetisjon og eksamensforberedelser.