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

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

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

<math>\frac { \partial \theta } { \partial t } = \nabla \cdot ( K ( \theta ) \nabla H )</math>

כאשר,

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

<math>H</math> הוא העומד ההידרולי הכולל (ביחידות אורך)

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

<math>t</math> הוא הזמן

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

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

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

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

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

<math>\frac { \partial ( \theta \rho ) } { \partial t } + \nabla \cdot ( \vec { q } \rho ) = 0</math>

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

<math>\frac { \partial \theta } { \partial t } + \nabla \cdot \vec { q } = 0</math>

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

<math>\frac { \partial \theta } { \partial t } = \nabla \cdot ( K ( \theta ) \nabla H )</math>

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

=== ניסוח תלוי h עבור בעיה חד־ממדית אנכית ===
<math>C ( h ) \frac { \partial h } { \partial t } = \frac { \partial } { \partial z } \left( K ( h ) \left( \frac { \partial h } { \partial z } + 1 \right) \right)</math>

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

=== ניסוח תלוי θ עבור בעיה חד־ממדית אנכית ===
<math>\frac { \partial \theta } { \partial t } = \frac { \partial } { \partial z } \left( D ( \theta ) \frac { \partial \theta } { \partial z } + K ( \theta ) \right)</math>

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

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

<math>\frac { \partial \theta } { \partial t } = \frac { 1 } { r } \frac { \partial } { \partial r } \left( r K \frac { \partial h } { \partial r } \right) + \frac { \partial } { \partial z } \left( K \frac { \partial h } { \partial z } + K \right)</math>

 <math>h</math> הוא העומד המטריצי (ביחידות אורך)

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

 <math>r</math> היא הקואורדינטה הרדיאלית (ביחידות אורך)

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

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

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

=== פתרון לזרימה אופקית ===
עבור ניסוח חד־ממדי אופקי תלוי <math>\theta</math> ובהנחה ש-<math>D ( \theta )</math> קבוע, נקבל את המשוואה :

<math>\frac { \partial \theta } { \partial t } = D \frac { \partial ^ { 2 } \theta } { \partial x ^ { 2 } }</math>

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

<math>\left[ \begin{array} { l } { \theta ( 0, t ) = \theta _ { w e t } } \\ { \theta ( \infty, t ) = \theta _ { d r y } } \\ { \theta ( x, 0 ) = \theta _ { d r y } } \end{array} \right.</math>

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

<math>\frac { \theta ( x, t ) - \theta _ { d r y } } { \theta _ { w e t } - \theta _ { d r y } } = \operatorname { erfc } \left[ \frac { x } { \sqrt { 4 D t } } \right]</math>

כאשר,

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

=== פתרון לזרימה אנכית ===
עבור ניסוח חד־ממדי אנכי תלוי <math>\theta</math>, ובהנחה ש-<math>D ( \theta )</math> קבוע, נקבל את המשוואה:

<math>\frac { \partial \theta } { \partial t } = \frac { \partial } { \partial z } \left( D \frac { \partial \theta } { \partial z } + K \right) = D \frac { \partial ^ { 2 } \theta } { \partial z ^ { 2 } } + \frac { \partial K } { \partial z }</math>

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

<math>\frac { \partial \theta } { \partial t } = D \frac { \partial ^ { 2 } \theta } { \partial z ^ { 2 } } + \frac { \partial K } { \underbrace {\partial \theta }_\overline {k}} \frac { \partial \theta } { \partial z } = D \frac { \partial ^ { 2 } \theta } { \partial z ^ { 2 } } + \overline { k } \frac { \partial \theta } { \partial z }</math>

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

<math>\frac { \theta ( z, t ) - \theta _ { d r y } } { \theta _ { w e t } - \theta _ { d r y } } = \frac { 1 } { 2 } \left\{ \operatorname { erfc } \left( \frac { x - \overline { k } t } { \sqrt { 4 D t } } \right) + e ^ { \frac { \overline { k } z } { D } } \operatorname { erfc } \left( \frac { x + \overline { k } t } { \sqrt { 4 D t } } \right) \right\}</math>

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

<math>\frac { \partial } { \partial x } \left( K ( h ) \frac { \partial h } { \partial x } \right) = 0</math>.

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

<math>q _ { x } = - K ( h ) \frac { \partial h } { \partial x }</math>

כעת, נשתמש במודל גרדנר לתיאור <math>K (h)</math>:

<math>K ( h ) = K _ { s } e ^ { \alpha h }</math>

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

<math>q _ { x } = - K _ { s } e ^ { \alpha h } \frac { \partial h } { \partial x }</math>

נפתור על ידי הפרדת משתנים ואינטגרציה, בהינתן שהעומד המטריצי ב <math>x = 0</math>   הוא <math>h _ { 0 }</math>:

<math>\int _ { h _ { 0 } } ^ { h ( x ) } e ^ { \alpha h } \partial h = - \frac { q _ { x } } { K _ { s } } \int _ { 0 } ^ { x } \partial x</math>

<math>e ^ { \alpha h ( x ) } = e ^ { \alpha h _ { 0 } } - \frac { \alpha q _ { x } } { K _ { s } } x</math>

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

<math>h ( x ) = \frac { 1 } { \alpha } \ln \left( e ^ { \alpha h _ { 0 } } - \frac { \alpha q _ { x } } { K _ { s } } x \right)</math>

==הערות שוליים==

קטגוריה:משוואות דיפרנציאליות