משוואת ריצ'רדס

משוואת ריצ'רדס (באנגלית: Richards' Equation) היא משוואה דיפרנציאלית חלקית, לא ליניארית, המתארת זרימה של מים בסביבה נקבובית לא רוויה (למשל, בתוך קרקע). פיתוח המשוואה מיוחס ללורנצו א. ריצ'רדסתבנית:הערה על אף שפיתוח דומה הוצג תשע שנים מוקדם יותר, בשנת 1922 בספרו של לואיס פ. ריצ'רדסון.תבנית:הערה ניסוח המשוואה מבוסס על משוואת מאזן המסה ועל חוק השטף של דרסי. פתרון שדה זרימת המים על ידי משוואת ריצ'רדס הוא בעל חשיבות בתחומים כמו הזנת הצמח, השקיה והסעת מומסים.תבנית:הערה לא ניתן לעסוק בחיזוי וניתוח של תנועת מומסים בקרקע ללא פתרון של משוואת ריצ'רדס ולכן המעקב אחר תנועה של נוטריינטים, חומרי הדברה ואירועי זיהום מחייב את פתרונה.

ניסוח וקטורי של משוואת ריצ'רדס הוא:

$\frac{\partial \theta }{\partial t}=\mathrm{\nabla }\cdot (K(\theta )\mathrm{\nabla }H)$

כאשר,

$K$ הוא מקדם המוליכות ההידרולית (ביחידות של אורך לזמן)

$H$ הוא העומד ההידרולי הכולל (ביחידות אורך)

$\theta$ היא תכולת הרטיבות (ביחידות של נפח מים לנפח כולל)

$t$ הוא הזמן

תיאור המשוואה

משוואת ריצ'רדס כוללת שני נעלמים, העומד הכולל, $H$, ותכולת הרטיבות, $\theta$. לאור זאת, כדי לפתור את משוואת ריצ'רדס יש להוסיף משוואה המכונה עקום התאחיזה המתארת את התלות בין תכולת הרטיבות $\theta$ והעומד המטריצי ($h$). בנוסף, היות שמקדם המוליכות $K$ תלוי בתכולת הרטיבות, פתרון של משוואת ריצ'רדס מחייב שימוש במשוואה המתארת את מקדם המוליכות כתלות בתכולת הרטיבות. עם זאת, השימוש בעקום התאחיזה לצורך תיאור הקשר בין שני המשתנים התלויים במשוואת ריצ'רדס מתאים כל עוד התלות בין שני המשתנים הגיעה למצב שיווי המשקל שמתואר על ידי עקום התאחיזה ושהעקום לא סובל מתופעת ההיסטרזה.תבנית:הערה

מפתרון מלא של משוואת ריצ'רדס ניתן לקבל את השדה התלת־ממדי של תכולת הרטיבות בכל רגע בזמן. לאור הקשר בין תכולת הרטיבות לעומד המטריצי, ניתן גם לקבל מהפתרון את שדה העומד המטריצי שהוא אחד ממרכיבי העומד הכולל. היות שהמשוואה מתאימה לתנועה האיטית שמתרחשת בקרקע (מספר ריינולדס קטן מאחד), מרכיבי העומד הכולל כוללים את עומד הרום והעומד המטריצי בלבד, בעוד תרומת העומד המהירותי זניחה.

המשוואה בצורתה המקובלת לא כוללת אברי בור ומקור שנובעים, למשל, מצריכת שורשי צמחים. בנוסף, היא אינה כוללת מנגנוני זרימה שנובעים מהבדלים בצפיפות המים (למשל, באזור הפן ביני) ומניחה שהמטריצה המוצקה שמרכיבה את הסביבה הנקבובית לא תופחת ולא מתכווצת. לשם תיאור מנגנונים אלה ניתן להוסיף למשוואה איברים מתאימים שיוגדרו עבור כל אחד מהמנגנונים הללו.תבנית:הערה

תיאור של שדה הזרימה בקרקע רוויה הוא מקרה פרטי של משוואת ריצ'רדס. היות שתכולת הרטיבות בקרקע רוויה שווה לנקבוביות, משוואת ריצ'רדס מקבלת את הצורה של משוואה המתאימה לפתרון של תנועה של מי תהום, $\mathrm{\nabla }(K\mathrm{\nabla }H)=0$. במקרים בהם מקדם המוליכות ברוויה אחיד במרחב, משוואת ריצ'רדס מקבלת את הצורה של משוואת לפלס, ${\mathrm{\nabla }}^{2}H=0$.

פיתוח משוואת ריצ'רדס

ניסוח וקטורי של משוואת מאזן המסה הוא:

$\frac{\partial (\theta \rho )}{\partial t}+\mathrm{\nabla }\cdot (\overrightarrow{q}\rho )=0$

כאשר $\rho$ היא צפיפות המים ו-$\overrightarrow{q}$ מתאר את שטף המים. משוואת ריצ'רדס מתאימה לזרימה בלתי דחיסה ($\theta \frac{\partial \rho }{\partial t}+\overrightarrow{q}\cdot \mathrm{\nabla }\rho =0$) ולכן, ניתן לכתוב :

$\frac{\partial \theta }{\partial t}+\mathrm{\nabla }\cdot \overrightarrow{q}=0$

מהצבה של חוק דרסי $(\overrightarrow{q}=-K\mathrm{\nabla }H)$ לתיאור השטף והעברת אגפים תתקבל משוואת ריצ'רדס כדלקמן:

$\frac{\partial \theta }{\partial t}=\mathrm{\nabla }\cdot (K(\theta )\mathrm{\nabla }H)$

ניסוח תלוי θ ותלוי h

הניסוח של משוואת ריצ'רדס הכולל גם את תכולת הרטיבות וגם את העומד המטריצי ($\theta$,$h$) נקרא הניסוח המעורב (Mixed form). בספרות ניתן למצוא שני ניסוחים נוספים – ניסוח שבו המשתנה התלוי היחיד הוא $h$ וניסוח שבו המשתנה התלוי הוא $\theta$. המעבר מהניסוח המעורב לשני הניסוחים מתקבל על ידי שימוש בכלל השרשרת.

ניסוח תלוי h עבור בעיה חד־ממדית אנכית

$C(h)\frac{\partial h}{\partial t}=\frac{\partial }{\partial z}(K(h)(\frac{\partial h}{\partial z}+1))$

כאשר הפונקציה $C(h)$ היא הנגזרת של עקום התאחיזה $C(h)=\partial \theta /\partial h$ ו-$z$ חיובי כלפי מעלה.

ניסוח תלוי θ עבור בעיה חד־ממדית אנכית

$\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial z}(D(\theta )\frac{\partial \theta }{\partial z}+K(\theta ))$

כאשר $D(\theta )$ מכונה מקדם הדיפוזיביות ושווה למכפלה של מקדם המוליכות ושיפוע עקום התאחיזה, $D(\theta )=K(\theta )\partial h/\partial \theta$. הבחירה בתיאור של מקדם המוליכות כפונקציה של תכולת הרטיבות עדיפה על תיאורו כפונקציה של העומד המטריצי ממספר סיבות. חלק מהסיבות קשורות ליעילות הפתרון הנומרי של המשוואה וחלקן מבוססות על ההלימה בין התיאור המנגנוני של השינוי במוליכות ההידרולית בקרקע המתרחש עם שינוי בתכולת הרטיבות לבין הביטוי המתמטי למוליכות ההידרולית כתלוי בתכולת הרטיבות $K(\theta )$. להלן הסבר פשוט לצידוק המנגנוני : עם עליה בתכולת הרטיבות, חלק גדול יותר מנקבובי הקרקע מלאים במים כך שמסלולי תנועת המים בתווך הנקבובי קצרים יותר, המרחק של חלקיקי הזורם מהדופן של חלקיקי הקרקע גדול יותר ועל כן מקדם המוליכות ההידרולית גדל.

ניסוח המשוואה בקואורדינטות שונות

כתלות במרחב הבעיה, פתרון במערכת קואורדינטות שונה מהמערכת הקרטזית עשוי להיות דרוש. לדוגמה, ניתן לתאר את התנועה התלת־ממדית של מים בתוך בקבוק מלא חול תוך כדי יציאתם דרך פתח הבקבוק על ידי ההנחה שהזרימה מתרחשת תוך שמירה על סימטריה צירית (בעיה אקסי-סימטרית). הניסוח של משוואת ריצ'רדס המתאים לדוגמת הבקבוק יהיה ניסוח בקואורדינטות גליליות עם סימטריה ציריתתבנית:הערה:

$\frac{\partial \theta }{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(rK\frac{\partial h}{\partial r}\right)+\frac{\partial }{\partial z}(K\frac{\partial h}{\partial z}+K)$

$h$ הוא העומד המטריצי (ביחידות אורך)

 $z$ היא הקואורדינטה האנכית (ביחידות אורך, חיובי כלפי מעלה)

$r$ היא הקואורדינטה הרדיאלית (ביחידות אורך)

דוגמאות לפתרונות של משוואת ריצ'רדס החד־ממדית

מאחר שמשוואת ריצ'רדס vht משוואה לא ליניארית, והיות ויש למצוא פונקציות שמתארות את מקדם המוליכות ועקום התאחיזה לצורך פתרון מדויק, פתרון אנליטי כללי לא קיים. דוגמאות לכלים שפותחו לשם פתרון נומרי של משוואת ריצ'רדס כוללים את התוכנות HYDRUS,תבנית:הערה TOUGH,תבנית:הערה SWAP,תבנית:הערה STOMP,תבנית:הערה MODFLOW-SURFACT,תבנית:הערה VS2DI,תבנית:הערה HYDROBIOGEOCHEM,תבנית:הערה FEHM,תבנית:הערה UNSATCHEM,תבנית:הערהתבנית:הערהMACRO

כלים אלה משתמשים, על פי רוב, במשוואות מקובלות לתיאור עקום התאחיזה כמו משוואת ון-גנוכטןתבנית:הערה ופונקציית מקדם המוליכות כמו משוואת מועלם,תבנית:הערה תוך שימוש בפרמטרים שמתאימים לקרקעות ולתצורות הסלע בהם מתרחשת הזרימה. הפתרונות האנליטיים שפותחו עבור משוואת ריצ'רדס מוגבלים למקרים פרטיים בלבד, כמו פתרונות עבור זרימה תמידית חד־ממדית ומפותחת היטב (fully developed flow) שמתארים את תנועת המים בתווך בלתי רווי כתוצאה, למשל, מהתאדות קבועה ממפלס מי התהום. דוגמה אחרת כוללת פתרונות עבור זרימה חד־ממדית אופקית או אנכית תוך הנחה שהמקדמים כמו מקדם המוליכות ומקדם הדיפוזיביות קבועים.

להלן שתי דוגמאות של פתרונות כלליים למשוואת ריצ'רדס המבוססים על קירוב למשוואת החום. פתרונות אלה אינם מיוחדים לתיאור זרימה בתווך נקבובי ונעשה בהם שימוש גם עבור בעיות מעבר חום, זרימה צמיגה, דיפוזיה מולקולרית ועוד. בהמשך, מובאת דוגמה לשימוש במשוואת ריצ'רדס לתיאור זרימה אופקית חד־ממדית בקרקע במצב תמידי.

פתרון לזרימה אופקית

עבור ניסוח חד־ממדי אופקי תלוי $\theta$ ובהנחה ש-$D(\theta )$ קבוע, נקבל את המשוואה :

$\frac{\partial \theta }{\partial t}=D\frac{{\partial }^2\theta }{\partial x^2}$

עבור תחום חצי אינסופי, תנאי התחלה ותנאי שפה מסוג דיריכלה כדלהלן,

$[\begin{array}{c}\theta (0,t)={\theta }_{wet}\\ \theta (\mathrm{\infty },t)={\theta }_{dry}\\ \theta (x,0)={\theta }_{dry}\end{array}$

הפתרון של המשוואה הוא,

$\frac{\theta (x,t)-{\theta }_{dry}}{{\theta }_{wet}-{\theta }_{dry}}=\mathrm{erfc}\left[\frac{x}{\sqrt{4Dt}}\right]$

כאשר,

${\theta }_{\text{ wet }}$ ו-${\theta }_{\text{ dry }}$ הן רטיבות הקרקע ברוויה ובתכולת הרטיבות השאריתית ו-erfc היא פונקציית השגיאה המשלימה.

פתרון לזרימה אנכית

עבור ניסוח חד־ממדי אנכי תלוי $\theta$, ובהנחה ש-$D(\theta )$ קבוע, נקבל את המשוואה:

$\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial z}(D\frac{\partial \theta }{\partial z}+K)=D\frac{{\partial }^2\theta }{\partial z^2}+\frac{\partial K}{\partial z}$

על ידי שימוש בכלל השרשרת, ניתן לכתוב :

$\frac{\partial \theta }{\partial t}=D\frac{{\partial }^2\theta }{\partial z^2}+\frac{\partial K}{\underset{\bar{k}}{\underbrace{\partial \theta }}}\frac{\partial \theta }{\partial z}=D\frac{{\partial }^2\theta }{\partial z^2}+\bar{k}\frac{\partial \theta }{\partial z}$

עבור בעיה חצי אינסופית, בצירוף תנאי שפה דומים לדוגמה עבור הזרימה האופקית ובהנחה כי $D$ו-תבנית:כ $\bar{k}$ קבועים, הפתרון יינתן על ידי :

$\frac{\theta (z,t)-{\theta }_{dry}}{{\theta }_{wet}-{\theta }_{dry}}=\frac{1}{2}\{\mathrm{erfc}\left(\frac{x-\bar{k}t}{\sqrt{4Dt}}\right)+e^{\frac{\bar{k}z}{D}}\mathrm{erfc}\left(\frac{x+\bar{k}t}{\sqrt{4Dt}}\right)\}$

פתרון לזרימה אופקית, תמידית בקרקעתבנית:הערה

משוואת ריצ'רדס עבור זרימה תמידית $\mathrm{\nabla }\cdot (K(\theta )\mathrm{\nabla }H)=0$ עבור בעיה חד־ממדית אופקית, כאשר מקדם המוליכות מבוטא כפונקציה של העומד המטריצי היא:

$\frac{\partial }{\partial x}(K(h)\frac{\partial h}{\partial x})=0$.

ממשוואה זו נובע כי הביטוי $-K(h)\frac{\partial h}{\partial x}$ הוא קבוע וגודלו שווה לשטף בכיוון ציר ה-$x$ $\left(q_x\right)$:

$q_x=-K(h)\frac{\partial h}{\partial x}$

כעת, נשתמש במודל גרדנרתבנית:הערה לתיאור $K(h)$:

$K(h)=K_s{e}^{\alpha h}$

מהצבת הביטוי ל-$K(h)$במשוואה, נקבל:

$q_x=-K_s{e}^{\alpha h}\frac{\partial h}{\partial x}$

נפתור על ידי הפרדת משתנים ואינטגרציה, בהינתן שהעומד המטריצי ב $x=0$ הוא $h_0$:

$\underset{h_0}{\overset{h(x)}{\int }}{e}^{\alpha h}\partial h=-\frac{q_x}{K_s}\underset{0}{\overset{x}{\int }}\partial x$

$e^{\alpha h(x)}=e^{\alpha h_0}-\frac{\alpha q_x}{K_s}x$

נקבל ביטוי לעומד המטריצי כפונקציה של המיקום:

$h(x)=\frac{1}{\alpha }\ln (e^{\alpha h_0}-\frac{\alpha q_x}{K_s}x)$