From 7be5558d0ab7719c7f86724227b3bedfff8b5cb9 Mon Sep 17 00:00:00 2001 From: Willem van Verseveld Date: Tue, 12 Mar 2024 09:32:31 +0100 Subject: [PATCH 1/2] Support different vertical saturated hydraulic conductivity profiles in `SBM` (#340) * Ksat soil profiles initialization Exponential and exponential-constant profiles, based on `z_exp` (depth for which exponential decline of ksat is valid) and profile setting in TOML file (default is `exponential`). * Ksat profile vertical flow Co-authored-by: Raul Mendoza * Ksat profiles for lateral subsurface flow Exponential and exponential-constant profile. Co-authored-by: Raul Mendoza * Additonal ksat profiles vertical flow `layered` and `layered_exponential` Co-authored-by: Raul Mendoza * Additonal ksat profiles lateral subsurface flow Initialization of `layered` and `layered_exponential` profiles for lateral subsurface flow and computation of equivalent horizontal conductivity for layered profile. Co-authored-by: Raul Mendoza * Kinematic wave ssf and ksat profiles Add kinematic wave ssf for ksat profile `layered` and `layered_exponential`. Co-authored-by: Raul Mendoza * Update capflux and deeptransfer for ksatprofiles Computation of ksat at zi (capflux) and deepksat at bottom soil (deeptransfer). * Updates and fixes related to `ksat_profile` * Add tests and couple of fixes Fixes are related to ksat_profile (vertical and horizontal hydraulic conductivity). * Number of layers ksat_profile For `layered` profile use `nlayers` for number of layers with `kv` values, and for `layered` and `layered_exponential` profile use `nlayers` to compute `ssfmax`. * Update docs `SBM` * Update docs lateral subsurface flow * Update changelog And fix small typo docs. * Refactor initialization of lateral subsurface flow based on `ksat_profile` Make use of smaller subfunctions. * Add ksat_profile options to docs lateral subsurface flow * Add variable `z_layered` for `ksat_profile` For the layered_exponential profile `z_exp` was used, this is the depth from the soil surface for which an exponential decline is valid. For the layered_exponential profile the use of `z_exp` is confusing and `z_layered` (depth from soil surface for which a layered profile is valid) is a better term. * Add `z_layered` to params_vertical.md * Add `ksat_profile` conceptual figure to docs --------- Co-authored-by: Raul Mendoza --- docs/src/changelog.md | 3 + docs/src/images/sbm_ksat_profiles.png | Bin 0 -> 71047 bytes docs/src/model_docs/lateral/kinwave.md | 64 +++++++--- docs/src/model_docs/params_lateral.md | 31 +++-- docs/src/model_docs/params_vertical.md | 22 +++- docs/src/model_docs/vertical/sbm.md | 48 ++++++- src/flow.jl | 57 ++++++--- src/horizontal_process.jl | 137 ++++++++++++++++++-- src/sbm.jl | 102 +++++++++++++-- src/sbm_model.jl | 27 +++- src/utils.jl | 170 +++++++++++++++++++++++++ test/bmi.jl | 6 +- test/horizontal_process.jl | 2 + test/run_sbm.jl | 106 +++++++++++++++ 14 files changed, 699 insertions(+), 76 deletions(-) create mode 100644 docs/src/images/sbm_ksat_profiles.png diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 0b1f22c54..318667d20 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -32,6 +32,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Checks to see if all states are covered in the .toml file. If not all states are covered, an error is thrown. If there are more states specified than required, these states are ignored (with a warning in the logging), and the simulation will continue. +- Support different vertical hydraulic conductivity profiles for the `SBM` concept. See also + the following sections: [The SBM soil water accounting scheme](@ref) and [Subsurface flow + routing](@ref) for a short description. ## v0.7.3 - 2024-01-12 diff --git a/docs/src/images/sbm_ksat_profiles.png b/docs/src/images/sbm_ksat_profiles.png new file mode 100644 index 0000000000000000000000000000000000000000..6e6bbeda2e4510bd77fbd2da497dea18764ac06b GIT binary patch literal 71047 zcmd?RgHxD*H!TjxZSKRp>v8 zr>SBtFxXzI)B|zVr+V`{RI#d4aqHgRRtq@Ku^73#6iVWXH8l1|c%vTO@_XO(KsHt> zSoR~YhYtLe?Ef$p%;f1 z`JZ<6ug^^X`3N&OeC7M0|NbK2?qT=;e)c)@#^CgSf5A<8W}f1|zj*k!FXF$SYmPto zbnL&spnvxNr%Cwm%Y-tJnDYx-)#|~h4$IFt2vCF z`0LB*tdx{G#FLMLj??A~gGI>v$A0cVAB4aa`tF}bq$u9_$s{Q3}o+U+qy*ID?~r z8N|-&@>&lTTGX^1c9o^8e%5#Sp^VJi)yszvSXW+-yK;$UpJw*g$?c_-cMQgzhS=Sh zL>5_wl}8jBH=dd~Je(Ox>&$D`?zA5G8Wql=JF-N~TUJ;1vWY0wX z+mUY$>>9b3UBVa(u(KV7Yq#+y?j1<#(fdF~xzlWhux*Si96BYPj5G_Q?n}Mxal&zV zR7C6H5+#}Nn+|Ifu$m2H^eVx1$fme7IjOg=PiO&4x}}CPU!eLunQ>J}x*G~d?d35? zPMj}E{_MG4KIVKUnZJiOmQ4%>Gnl{iBeqJm)2!s#4?e5DSjD`?7bta%-ooO(8r6)G zd&mwC4-Y{)#NiGP8wVlM+!sq58ECyuNNEJGX{$ZhJs1hY8j$3RxZYFXG^Y&z>+wB= z0)iL{o%!gA{aJjGgwDd0C`QGv*b;Dz+2dD=Vo)hCxwCwVnomEm`Mf~wKN%}Wo^P;V zrYTaN&|7R2BshSz>}4BqU#e^uO#+sUJ92jcaV)Fv{-5a$JcN3`POQ!!)lswaMTCWg zp|+PB`5dQZXZrqzw&+!tw{D)nXG&M6;WW7z(c{EO=VMS`0GqKb2%%;u1k9G&H~f6i|`x6mvD1KyR`eWXg^Gzg>h~;KYWmM z3qD{yr0(ik%=C1eVtv+OY-@CV_zVnIPg>t<4$fzD9E_rsnj&6QCvc zQee@2-PK8FNeZiYBgM`WuU&@EEsx#No95a1au*D^!@d>?%hTc~*mamhmLAHvzdwob zFOaXGcO)hF@Fp>JUrff>d9B7Q&X!0rd-gDbv%Yp;WT&4lZ6zV|aEaFr3W!iE0(Dsp z2L1%K{nKkFOIug%%M(CB)6^O`iial_SrBXAjd62x4(p;7OJ3{ckhMsYjk*%lIBA`l znu;tirBQnZ{Fww<%mR7U_QG#!O}s$KmZW>*Il7~>u-Lwlef#LcOIGe_f*4gX`B*;= z*_Hz}kIGt+E1CInE9r0A+0vPKcf!VGkT9Is>L=E{+Y-R+_VWw}vT!zqeGe%8*B);< zL*LIkS+G&pG0(*tP1(mC-|nI_wXmz*S}rk4*D zjR=IhqSfZ(@lkNv0+0D@00r znw_2Cvw8uPN5gTZ-5Kc??eA+EXNXNt@EXkg%Nz|`JI%Ylr}XaRC5K+nd-|D?g$DQq zV8J}^!5TM4*Q;g%Up2<6;%!O|vZ<6oY* zP?r2GEHuX4-4Wa+if;`u)SA#2yvq!rfb#Z!diI#|`*x z7hjd4tM;Ty_O1Hlr7?R+sp~1mwL)KraZ4VJj5J;>-k7TO7?+DWCZ{2psiq#!tgYSY zr)`*f^>hm(E}M}qM$3)UKGK=(i5hyg56$-XGT@oKWgMC=!ce}4m|!;OKve(1;9q*S z+RiyahO@Mek!Ifa1Sc2~GI{h-^2OANS~3KE1wsyv|8iQ9*-}}pNdM2KkT?hL=^{eT zUytbS{+_+TJXuY7N%~rLgGMJ^?qyeN+g0`g`46$IwBderOED%eiGXdNl>mm(=fey* zl8V)KVc_K5;=3~`zw;;Zs?yyp7@fUz`-LHZ=>gHWr{oXDdc=7GTk+QPNRir6me+` zWP~GLTkKe%53dP5{Wb2BG@-ZCe{Kv|Awtce(oqOy;V9vSLAbA(Hv}*& z0^d%LAz{wQnj}JI#T|}kKA-mIfbY1CXv1K?{)xPE(No)^;#chnAh1JN8%N2mSV2cJ zc9Pud)8XLY0KZ2d5Q{e~GH{~%%X`T6dwGgCobCa?p+>JQ;^Hd|f?{T60PEFrkf{KJ3v6AV zH~xp^KGvg8h6LOG{D{}q0%|_5akw%hN4H!(O9SJw4#t*$OetI$wkLb{>14w@ zxauaa8cc4(tUFDvy-G_Xk_biNp!wwaK_RJGROTpv=8+38FX(jwxF44_D_72$n+CSx z{y?c?I?zsBHvD9P*Tyt|_1p`O*1-wkvCZ~lHYKur2>_c2*i`)P4;Da%mEZ%txNJRe zPG(D^m5U!jApDZqbKm?xDTkoBgT`GO>5sT{$GS6Dw_G&iusM27HeN)lcodaIdk`6O zul>J*WQ_X(-X(WpaQGi^;Qwg>cdy?0UI28??En0_;c8Iac%WqHsA30NEh`pCCe#AB zV(!>HxdVf>zua1?9%D4DFj}C12}EeR8 z{uBKX@a)ZpjZA%KVLpz2ML?8imrjTcE(8X9c{SdIMEln_gPms|eARV+&++$9zy{O< zeLH)FN2bW;Oy!1cmfpY^8XF&Jv|JbTDCn^k_H4EDLV(WB6rTz8#@S;4E0UzM=};5l zr}$qHhZnCepmBsQKik^FZj$5tQ2*2pik` zJ~gT|p!_TxRESSV$mtr}%T?y8>i^ic7-{6U)w+{cS9XPe1cNu5!lezA{s~4-UP`!5 z32R`a9n7?OH#@bCB>qZc8(EKej`Fg5TboxR!W496`Ef@2HQ@}2850X)7z-)p5?g2j zlYIhMWPoehpuU>pgG=j$v%>h(_$6;|O}opMu6<|HZrTbsR5u;cc$_{PJts&Sbf!ak zN@RBVR{cB(2Wq%GfktuL66$U;+ML3-U5!I6D8{VyI}6+Jf_3OB*k<`>q5=eE|7d_};DYX<87XUu8U;s=GqyIG^ceW}HP zF6gV5H7sl3-Ae4?8+jM}2H2WU&iY#nMC*nxs^X*+mP>O%AOm|Up#x0;v zRmNzlJea`V^TMbQJF;1f9Aswb(r_!RTG6`kxfT+c@;;E2VWyGKrU2^oj(F@Lt>=~Mbd;EX32|^w_%XXq`;+D*EXpjASk4=A-95Fco0U59VE(}q z4%!$-%l@@y^-BG*Mhdl$dj_9RGqVEDH!%Fud(1k>I{_TJGQwVh^-zi|4YDydb=S(f@I5P9D zew!1FB5IOohl&O%p#d$e<6)sW!kw4C)cT>YAD(-ILa!Y?s%Arok?=OmObU_-+NB&bq7L(;1mgYg%*_W}( zx#rV?NrsRV$o`KcUG*(Bsw!id*HN;gL`X)X&QU*|B0e3awGZ!lE-Kf4+QL+;?7UoC zU#Rc;H_Xw|*Yx6C`(EQnjOEB<1z4W@vZr7_IYn1A^*fKH)9qoRO5OXKq zH0O+lgzckZT_rR|mtP$0P|Cp6oU|0=dB`o9D#4Nx*3PIsNb<_s-7LLInq%i1bV#{B zTBmTuF_`0Gap>7!z%6f*hAzQz-4W+4)fwR9K;iuc%Z*7RG4Dvbj*`SWD`&Lzadvac zTD>3uL?pCW%ZT|#p3WAfwpBOuRZiSr*)lKQM@1vD0x>UKDwZm!hUXC@*Gfd&a=gY*RIh71@DY5T1=mX zf^&Mee~aq&Tg~uW%G>hX8~FY(Jg?snF!$}Azn2H&j&|rf2mAFQ&)PykE*&h8ZD*(< z8S>_Xs0nIP2z+2$FKAk!xod1#KwQ5Jt;#FY6?#5p>O2^){lzo=!7=n+!H+k)j#fFN z=tV^VTh{TTFL;$bkfOr+-o@6wyViOi+Sz4O?;S`F&vQt6t||eqcv@csy|lah20L%3 zrmfIw@xk!$ww@g4NkNWVRW3F!If40eq6EQOVuG%u5TSA=(0pMUS68Bo;E@rf9hBO> z^CRxaOq7rP%c_zL)tL+Ff`-idtD=`7%gM!-e6{}b-@L1kSzBJ7id+>#AHx?*99P=% zQxeCrSYQ4#;v(UK#9U_j9`roOMw#L*fhn~5l(?l57 zil&9l6xQDMbe4)zNi}PKA$!C#8-6gFjR!R~Y@<`wzEO47&BZl@X?y3N9V7Qjp<0I- ze^Q1?3pS-`^PinLmQ9j=3w%z2lU30_ff?-Oc!~M2x`}bnuCaE`?)s}Q0LsMF@)R_` zzv8#0yv@{hL#q42g|JyWyXZu>5MH{!p&3HWAa=Z45mxiMYml+qu{?$t>7Cx`U;L@(*{ zTS@1PyvEYt(b>QTwMe?omA5N|FcuS1mzqiR`@)jPY@ zp<#i4rha5cZ|0bkhCpt>xrE%94_U$(siP>)8&V3QAop3L)_wLHsrV9+>rxhBG6Egv zV_$=JW+OOB-c_P_QM5qY>0iwqJ{EpTUL%) zmEn9ksel>aV02}gof>;nqxo4?l3-%$U|kM2BRq%*&HEzJnhM*G9gAGlidtrf0yZLY ziH+He03UQ@b{61P&dPDgtJFw-1p7{ID2JoxVfMja-XQ=6Vy>`P?O3o-H#Td{Y&9MV zHW>h(HmMwwkXWsUOlnDg%v#?O1VRdXott@j-}MI4Bx$Yu5^OT%0*2yTv`yO#C~L1* zIF7XH?6lF)k!ov=;GL6Xl_3uW{~>!5HhzdhmuF`y{Nm1%*(sJo2u7aN4NrOYleDv} z?Wbm2G{WZWY-1DI(@fcg+m^^%2${L6_4zZ1m3u{4Mr##~I+sCP*1UM`q*-P?K-XXyZ3qzaC9&bdBz1@ME7{ZJzuf^u@vCg$Fpn-1-YY%|D0U z+^p~eOb(26T#{n}1SUi~KFj_Or|Qqzv9zJCTr55Z0+XgwxC#H0K(jSKw9#yJ9zao8 z`rQAG+U?i!)?nNhq4_qHBy$)diVDQ{HeP>0;ocpUFoL4(S_heM8~jKAlR0G-n(M>d z;u^ATDFiHUO*O%L@J^y+x+M;BMMMY z)u}s#*<}qpVZ)ia&53j@q`@xgwy5_>cPGF%$gO7J zi&o6b$l1!}UjQ{x<0-xiumm5TDE-#rpVvHk&k zS4uMa3u@&ZYNusYk`NO|w8Z-Xn8bFuAao<=E&DCvO9_@T_PqF#C&<7qUs3W`pb>ov zAZ{w)Oj=G~*0CWnLfXV4E7b?43!5HMJ?*6}2#sLsgmbp{YlMf-6|BC!o%(mL=QSAY zuClSDDt^qe=QMy}d#D`w*B8G~AC)z4R2Q+u`%b8OdOf2$9x>OTwb@SEr?&y&0vS!? zn#E{&&%MoL((?308!wCm;!I2FlhGHR^L1Ez>F_F)P3f+YlAmd4O3VESVB9+;B?WS4 z)^Fp-?lARwTMMTliga)37Ws-|y1m$i`4YwAv=cG9!r}=PDK(j}uc(93L{T7AGdP!J zx%e%JUqcxFLl0G~)uzjnDxlFT!_xvm>MAeYm$bMjAfrBl_km3s_pb}4+hZElO!Xd9 zPG*!G1mp`f3DJ(Mjd|MmbRoccXGD&cFyJP^)$yE{Sz4aO%?*(B`4zerG{P)U4i_n)-~nnI^nlLyQA%My+Yh1?;@;NMedyE~ zZK&W&zXYn|u=L(S>38=g)1=D4N#D9B>Q;gV;SND*H$`S?STy<2Y52K*yh|i$*Kp1_3pfm^l5!Ah4KHe=|Do5ABjKi?+-1Sc0fA}Rp zZ1Di#ypOqpVv-E3s@wZ;!6)bvG}PytGJ`7JnUHt$-uQ9`@>ipVGT$O|Vheq8Bg1p5 zOAUQ@anZNZjZPl{X7H}Q4T8sJxIY^}h3rQC91UIfYYsihktV`nJ8Vx>3MT+|!cR*f z?+7bI{o}pzjJ$(D=7Sia!ICV8uAoq6al3Otl7y9A^-U3`bee~**K$LBBqJ1>oN7h% zwUi;QqX{DUO0GD_xLXMmfh#3W)~OGC=Sb*3j7l+|sQTFIV;dcpUTzneq0qtM#0EPM zz=!7{v8Z)=H93e%lsj!Fgm%XCWhTQpt2>fE37(kwzzK9;OkXn_qpO|d7IG_&7OLs0 zly%QBwmAcGs>xM*!%8+IV{9$*H?N8SuRYN$wJMl{?Yo)NdEtDhxse_<-QND13|>G< zw5wm;oku{#QACQE$s)9WlLu^58a}*MNhajPJ!Dbh18>if!y^Lh`WE8)E8^g*xw(&Y zV%CdsTODe1e*J-&NmtQ#YtkRJy42Q0$hG(B>*-Eg7qoQ6{o+0t=O=lbyPoc<)dK^w zhf%wq%B38J)TL{(>|50>CtNlqHc7ayC!Q-49#iRCBa&e2uGLRq*ZQ0UDSj|lY37XZ z#Zhlm-g>$7k%!DFt)B2NGL?$rObAZx`P=RYcphp7x5vA6fV9_hXltJ?#OcS2eN3>+ z`u1)5(=BC>TPA>1kSS>CCJ^y#gbOF$@?<#44TR`XzL_xM* zjP~3eiXRvbj{cPstD$kyclRlB*h5-yM89}4El=U~@2JngW2^}3YU$fdm~K={Ve^fz z$1RC;Cdh&k9Y?q+TciNU4tHv_<`->_oZo42R3r6XA`11LSGqACW-`?738=X_b{RG& z{Pq>O*3(T^%@AEhZV&0&&mF1}JpXQpL~K9mxioMUW#XNMET5F7Aw7nB)!WGQw}1UM zHga3M#dId~2K%fiDEk{Pfnff!3?FGO99ijE_~qP|qi@Z~hZ6r=zNEKDSjw#OEXu<& zom6N5up~EbHo)ThIM(H-GQJ)n?rpoK#(PxCzm7)zTt`hUN49HYfl!ODt{i4{cEcbwUO(N$pzX){`72w;vg!nTz)K zzFaRcW?5|MI6hmjdcwJbhPLomg*9dE!Jakz%H)E19G118S-_n_l7S+0jr!$GP~J~jr4A8{te8B=wO*KbdEbaOKike;di(IO+SA7CgKB>zk-H* zOV?hBg|B)mhpG$PAG1yW76h0zT8{@GU$!~{T#;;6KxNU$%9cp}P3<3%sSu$h#;dWe zBPK3m`m-bmSP%iu%SEw1zfS83q4llB($d-^L12BwujqxAMy!4?x>{2)ySl7QM}=CM ztYhTWG+Y7FY#;6<$x%tL&Dvd#V%wR)8sknj#@=C=eyWbvf7{;Drx6~CU(0G%ea~Re z3@8Nq;cwG#c2|3a5QNukl_aln>mgu3Yrin5{Y;q+7d%3W!NV<^0(4T3x(Oekc4Nkb z!-r}ZZnYJ;^gGQp#$*v^Pe6f|c9M^VgA6XPjxUE>O~2x!s(y9uy$hi!EX)8KcS>5h z$^uZ+B_`&rS?!%Au6sttt6F2^TZuAdS8sPgY*AhM3>_QSbDslS6L0}I@bQy9SP8;L zRwDsqD{O|lL2SLURe?LH-7H0%<6WgDz5~<^kzjpfhiP32?1EHUW#1(r{a-F*5x=t+ z2`P?+X!1U#E4nE$8P0iArk@dDd@Q>(U(pS<(K_V5`+;$LL2mW*TuyQv>D$|0QZF7T z_C+Qh3cWr zQw^r;$%&u?hM=WrV_U~0Wn&n!f#lWzSsVN4)-a;NLPr0WOId|k==&7a8eR#KKkAnK zhu;I=_QZ1zeM?s}KG!G6+Y}V`G@a^kwOW|A1krjNPWYoVO{C8!4h5sRPaE6T@Xpor zh;tdYh3n59K8QX7RAL>D`gGf8HZR6la~^J zDdVoJs8-t}G+rHC@rvoTWuIODns1}_Gi^5h(sdj_iK@YJLmYs=`$WQK^x0%}D-`z# z85+|5wRX*Gvw~^F-Y~ujTd1QU*P?`c|L45_hM;I{Z65Xw)2NfGm{%r}kHau(5t(-N zVa?39E{0Ol^omq8?>}Ek2LYr>tkxXv3SZYI0-|U3ba_k^rHcg%4GnJ9LxlCovN3-8C!~Rc0ks8m{*~S<$E1bbniXLYC)vU?&q(+0)?q zn0_%4#|L3Fh?Pjd<1Q4olsGBBs&ZQelm)MX7js1KA(&4{J%@#>b9&X@_8H)>#DUdy zLtu?ShRYUtw^&a-n@^_scIAy;&x2!y!6BQwb2($FFPZCdFPA<>L%^~rWFz9%*@DoJ z^PL=yC!p9xuIY3si2p!h8fk+lPX}>oMqlpLIC*yt+Q06o?E!#?2SF|_l9Z-%BH#NIYDdR3KQsFP5h{YdMdge_2!;WK%2 z2qJ~N2+Bo$RXRC?*f$`!{m7#&nCZ3CRKH2GG(zxrKe?`UfKFgv?txs`Sxe+qQin|= zA*78Dcf65~`YD~CHuWX~!e_1=oc*y>42#yFciNZ6;aDc~?yE6@xKvkdz>!Q8e{|5E zp?gKBXP1ukv3F4K{C8v0LA_*6$@v0n8{zM1-d?%D@7o4j50&~>30qnydp+_lDv^a4 z?IR2n>!2D321i>3&0ZM76DHgEP}RBG;?@0Bgu@IzL-TdZk6OFd8+P%tSz-dv|XO5D@MD@mQ>M&D0gRX5p9Ai}?!Jgr$4MF0SoLd2b!AY5_z zMLNwWh=fy}XT*d|p+ZCk@iZfSS4S{yV_+rWz^uOZ{@4 z>C@+KybS>~gu!1+ea?bvC?@YBdFKVCm$LYrx{FtN>$9|{&qD?QgsZG#|DLPa8--|! zfJ@FcXFx#5hS+@QX>?f)eOzne*DV^G*@j!7HS7jv~^!MqnnIXI^#Z?GP4F0l5iL>B703y{RFxVN7Rh_jlZ< z;5F1a$K?VjUMmF57^Ce;Vfc`2XJf)fzugZvEz)4oRGB3Tq=%vOc}>kIqAwN1KE2Fi z7V2gK(}!Qk^;Ux%g(yGfk#}Ck1K@?Hh~)BP@(8g|IbML)$QcMUa-47?x^MwmIRE zdc;lppg>~7{1?@2P_VCP98^{C<2}GcA8gzo1N>;AZc$d{p!?D`O8akS4amTVNBaD- zbLE>o#-)x|y#HN`zivYLp}&4Yv`Qj|gSEf4>m)Sx@5W36{8B?`q=}WgN%^kLq|!p> zqEl{S&=8i=YcASi09lZt)BU!47^+=>XW0Su!7z7YMl%3OC6X7&O<|}_yQ)H%IpI%C zccc!Diq?bxAG!TH1I(5|>#QFVxu>4Y_A)zL|4+89T_2a4W%kU z1q4+6fMz!mRNJ#wA!2UI0VMGNBFdzV-u*FHr8V#_Si@o3@y+@*p&WN&G+zPqYIl^$mAyLuL`?~U@dsKg|;+0mD@?Q7&o>bf_TdW?Js zego$vAh#%kJoSdDyZ(G!&X`Re~dc@A#s|0yf<*4Jw zY?U?Yl>NuDfbA;$c-wOmeG9N{V_L@h=LP%a04N%Ro;g1xDJZ9rkl>e-i@?Rc0J10R-Bj&02^@miTN3A2| zn|sDaZ!0-oH!b6Z9AaM|V>ib2_Ns#}{^fl{r1j zyXp@ZL^ z{2Sa@aSyH<%dBE1!0A{8j2l!?z@@8(R?rFM5;OY&EY05^9fmA^Y$jg_!o9n{5WE1= zgJBcvvm}PnvK@aTIgbDAEWI0m1oR&QZ@(r}9-y!(^=L)DHS5_3&Y0LtjJ|#@Amjx# zK;0GWWQwrS+yFuj7+^lMh2IU1Sx}^PbIR;6p>$IM>E@P4&1*5WfBWQjN|xRMDorVh zu_~<0*A=hUK=5em1xfO^wti#Oh_HC#g%$=h|M+8oifo>iD^-Ia_FWv0TP0IEZY56Xki6|NlYLI5SC3pbXS@O{RZ9$XvxRO(^bj7DXbujy(L7v?%HgbHX3KvPZ2 z5l$Pc?#<9&vQ~6=)e#Mn#B#JHgo3Urj4ZLVgacb(a2r%V7}Rmop}g`-8#2F#OsunA zEC`Fi^(S!M=!0{yDn5L}jtA&@{G!km)Dy|dzUq;T#YZ9x<5JstAd9S$X=ZA1F-PQ} zzhr74QdfPch}rH%>P!93tmCWA`8ITq4?8vTnIJ{3_QwGm58(1zMm#E6yko&uJ>n8> z*)1h?W@h%piMML;F%y5($$$XnHg3G$kmRysZ3S}EU#O#{5e_nL3=~jR8j_o(MaLYN z6YdPs+8eK;WOp`LY#j4jsP5$$7Zb?)>QG?%En;|bxKl9tV&Wi;TUr1f zJF4D0oZP6F`%Pwvk3BE+y*2CA`YeTZ)=m&p3aGgs^hUq1QFy6F&?Gl)?3YWj7fJeE zEswS+M)E^sD#(B4XelQFzHd<$Ea4fW^(6f2L^C%(`7JB<{j=+LexlOI1cq}>gKKC5W; z(Q82Ts$@vIF?fyK;@9WzXm><99Eb|HI*3K>3)=$pI;p8C+aDm$%5{1Zt2QJYH?9=CGLn6-IIG3n=P?TU`4Q#~tRkNY4 z5W61y1x&A?vGB53UWpS_R0twG>4Hr1y&`ORTDQfNMJ<9=Jc0HvPEaNWLOG|ARuMOZ z>cQg%2gE{=?@+P+EE6?O(%!)(%J>A0B^5T4Vw9+Ndv_rXQuKU^Eu{gf}TZtD~b2 zWS29zdBdt1Z5r+sDJI~vLLKc{TH2^c?*+2k$$|7I^UdsYiL_y$454~Pl}S{&8Wq9g812WWX&)Wr4C{uwVBOF|_zMd#E-V;^E+=_|1@8#Lbqs zd^|z_9_BV6vd8bsyh2&VNh4V0fwOq%uaFR)lS9M9!rW} zAX?785x`RkFXOuaeE>>%?i~m!>>8x(gMk`23#A$$&)Ut$S=fn@il-L<+ghUtn0v0u z80+fZ?@Vg`9qdE-e0Iud>SzBb_*fZw0Nq4X#S`_TL9aEz>_t&9mvt77-UgoqNR-TkAVJ6e_kcCv2giz80w;s zQ9G&@o7!^SrzF*CKpL7G|<>MfR}dz8&6 z;y3kWy9rwIm4BCf(g&QOisdbEuPSq<@;to8I5%721d3R?Ne;F$i4Qm`3#%m*tT7sK zRKyKht7H5<1;PJL`#Cr{gGbkaRF1A`S3UMgFIMiLDS$3I_#{Osc1(jA)UToV3MJ3| zf{(puO$g1@p1F|e`#&8hHYs-SzUIUXa8Si)YvkBDnp)7EoyntQkfPyp{?^-d3kk@)PRkyn!A1ki$4zH(3kdGGICGo1@tA0pJP3Gq2S zpd+~bl)A~x85f!q7IDFf^55CKf@VYOBsmG~cO^kD2E@A{r!vPXx-$U8kDF|`=#j(u zkiYRFK8(Z3Wi&2*&_c8(M5bzUaRKix{p8BQ+5}j95in5$`8otpWv@N;5JwZwn|TRm zRzTACgmW?#@me6ztdkAW*^DliZTqXMm6OrapOg=!ZH7y(!eiF7Uq;{kccI;Y@2NSq zTavEj{WM*Q!uwW_aqg$Ro>ZYy=+ zT)IVl>dWZ&lK)QSG*r(5_^dSj3Uy5GCzc>M!@QONpA7JmH|jgYfMqc-=gTof5Fb8w zVVZ+v*1wShqB6r}PB?oxL<|491`*5I{{7J8FvKnwa5SopF)4MCz_i@NZ3uceaTPr$) zs;C4Ftc?|@fM@68Qnp3a{KkVI_KaC85=dhDa(k(jdeiw!s}4tW^a2#u{%*`npP39| zfqB7+6*G1(;yBYs!l{TqQ;tFP0U(*;YzCcLw(!}0nI)|Gm!Wt5^&tx)DwKL7PYrWH z{eYrZg5TyrXFj~(CKCPYLH7yNcCQfN)2>%50NVIT_cm$b<5^I|`WtGU_tPF!AIQ~T z$k32+1vgHNktjVrgf@R`v*SI;A&}d6D=;n%QIA^t08!AN<|!w=?0&8|1#A*6uGB+1 z{P$LYI*lB({LYgt9zK@`bidqwuA-cJ3#~m}^gaNB00Lsan|&IvvXgHF&#VvW5=l$C zH1b;zLja z0T9b&+o7?}6A8JcnAm8YqQV>Hx_ z|*J87udJqP!&+zS-1u`&vFh;v40%pWcVIyuYn~RYpR&3 z?&U2+f~yefHA_kWr1`o11C_{JpML#2DB!hvE$6OHXcd<8kI{?I2L_AS^mesqzIMB= z>OZm-1R@}V0RKAcoyI<6)RuWVr5EOkk_^iLE9tES3)%9s0=7dh%eo)q+Q%MDza_%7 zS;++G3I4n(*1SJjC~5>WhXPCHdei3vNn?NY6&nJshuaS1QjiKC|etMi$2 zJzM}xq_XGwbX>(WCZ}HbiBoYr5qfsq<{)Z>%}^$MJ!IVJ2r4OFkGxxq51m-YHJn&2 zGPlIKgNF?4!`_XcR){E2vVV>`x=RELGDzL|Z&yMYR^aeIyB>5IS{V-qZ>wKPAd`$9 zrqLoK&*@>HC?9+*4m`zO65aR7D@s#lK6g>ICUYAu$1dG%Hcf}V65$fRe{|ro611$J zj{s3qnXKo$a3!3YDljC5B3DKy(Ye3;G3^}e$p?=JrEE<`$k zBwIOjwB%M(S{vnM5MaLKC#=!GlsWeVX$Lko?c6H@z(8Ra$1y@sSHR2kK3_n&D1_V; zcBsLIF|OcWX_b~)E2}>NrGvklK)S%t=>Q6-mJNhTki*LS3~mVBU{B=IS4+zTq0&t# zR063rJ-2E#u`y&;`=LdiauM+XW7iVer8yys}Q~SM6yOQha<5Jba5;LWJG>9`ZLO3h;VEAs8*@YDnR2T(^&<5drmyozzdw(@8)wRslscL3jDD~ z(y!>4&o>`+othiH-_P`l?f!VP#)wz8tjWr?*8!CuUHY7yDqep~aNXdKYCxnm@ds;M zaYHacyhgnQ3-YO6xZ_Ev@1rX(nHO9nrNpX@!=bCK8|I&OYo0!=XEHd412ozXnO@YA zX2d-TOf6LKkctCIO9&YIsNAyhyD(<+&2G&4z2D|mB~xnphSQg9oqvH+z~4Y7lp6AF zb_ih?7et1WKo#<%R%I!5#f0h#O^A^tN9AuDY$`$VclDss%wZ2oQ3WHpAK#pLpA=fV z#Q*EO0uDocP0ukBQOhIhZV2UJY#0+k2>t^gi@A2N4Eji_<4fT315lc*%&TdNHE__8 zgM-9FBg6+xPu(BGPA>?|Lh@IrT5LYv6x3(97>&&=gG%02qTkyS!N8xwUMZ5RTrWDM zWJ>;1tq^Wv@lXma%3u9xlc8|+K>5hTRlr%|g3^wu#&bt2G*J`lR9|=#fR}r&fL*K^ z2Ni@M%_xCU!?KUBvYBc}zrP@z&Eg+-5IBPS&j36>^aB~U?y5oCnFp~JpZN-It9jqS zJo*exxJFU7sG2`Zh)~XdDUKiJ(DHgDSK4+{!S!Lxf+tp+L#P=j>B*bfH#NPH^0ngO@1&YlFf4PJp%Z%GmSH@SN3tdQo zu=9WT`hcTQ2k6q4*f&-^dDsk4XYtWNhF>=_{X|RG2m|QjIo_{<*u2F7E(4Vp*gdK3 ze%s}=5O8_(#?C=?KuyM7uop?#seC+$cS%gmf%`YW$I911D-j;N{LRrF z8LIhEKs>%SJG2gYMUt-cA85CKH~KjsWmh|Q%WM&q$)6cw=oJ8IZE|5M>+h@MWL94q z_@(MmEeyCv0=&eLZ6^3NkNU#P|LK7S>OmG<)Pnd!8pbn6Jebmv9YwrPoE_X=2?Nt% zU_YWH1FD}aw~|~y#l9P&3zP(cj1m6us{`9!GYcUiL%kg9?)IyU)JHu08#thvcbbBg zvs>HFz{#hp_g?OU?&2`!#lE-Kd!-LD>%DF-GZ)FW0a^?MN9{`_Fwl2whcl+9O$V@^wJyKL#b4h0 z$JGUT`I9qvOH;XfO_cjwrFvz4wc2!rwn7isI8~178)M|eyul4WmaJ3!YqO{A_%k!F zP8W`Y+ErYcR%tZy!3oZTAs=Qr=;xUY0Uk*i(cwI;7T zPX037Mrk|;1rzQ1j~uVG7mgua4C5WaHIW)ckb$;ScQzmTeCLQqQ||*eTYsp=$Feix zKb9obfnk2m1ig5Yu2U<|sGe$7Mq0YiYljqqCc(@m!q&Ax<(N7!nR#{Ukqn&)IaI;a`_x$YKwM<;D%tmhum{y+}72f-48sAh)2%Uu9fp*qjdx==j! z`)yNCDDPhH;dNViynft%8aI>-ZV}Po0vb%7h_F&X3DB%nfipliUnTKg?w8rpCdkXB z11>Y)M!TTpmeZ7~@9e>m1RR;<{9a%^aBBnUMSgPBmX|BO6o;aFs{cd-`iM7M2eo4^}PTc|gQW~=_^sDJeX zI1-DLlX(D&hLqakl-y`rUi>IEF0f!NF$an>!vO!d2}*6*bSdNAA_K* z{Cs70^=w52g&_9_@?lIWImZNsc7Mqy16RRTzvMh0k*q7i_>$11{iAKTL<6oeSch+{|GbZoPE|_aj$#bTX;gGkNlH0 zUhHwUcCz5Dlpu;A@PAPPqwyazM~#xdaUt3IH)}Nm4+<0-J&u5(-5?P|0Vy=6lm$GgHv|h&rmV{ zRw8d+KLQ+9NnGLq^_0$++@o;^=m8whhk_t_xHApf$ce8^)P5#G;i%@qhgjy{ABM&{ zoh%bVOO$4|!5g?ly0yC;V$>ja(wWhiH@~r3SUJzU$BTH~y$~8kx{e0tb1i^tW>21z zdeC%~jRU^sSZu!JWF9%q>;w#CNDP*}57Ts$?*4Br%rIN0npLkdC3_y&x2+R98tBPq)Qcv5- z?oR72?(_Li`=3LXCG4Xd5|span9#Q9!le7e9!GBB(T^ju?X`4gOFsRd$jhhgmTZ1sS7Q0^9DiM|?W62nwW|T!sP~bP=ZD=X!wbv; zc^o6(pfJC7JT_X2%nNP1<16ya? zP55|$oOGsPRxoup|5y?X9{;(2pxPpzZ)*?Ef*lxX9STREfOS^d4G&4Fqa{N%>+5 z#NT7{kIf&$YnP3FdbBl+&9nrX`a=ILDkfmaWb_y)O7fTWas`N=G?9aFHTj@LrL9#D zh$l8?2;l(UD+sA5ueTD(Gf;g65fU#W-gNiJ)7Ec!kIg|~rOXTu@B!u6;qG{eG1k9@ z`<}z@!w+N^WEHTODG$jOEJuE3>A3X2bwl{4qqkB_{g7#}KNRfgLJPX5Yxy>2!8iA8 zYa&IFIq|)?`TL%XBDeQ{Zpl_9`80Gw$0|7b|6l#W`~7f@=LO*@#~^A6U9>zt@702JOY&)Cqff7Tc^ps-eH16kWK4FHBm#65$ zXJZ}SOh9pc$H5Y`QOaN-!k+XR#83|Y&$**OQqHq84`5WF_A6v@oXWkH$`b{sK8m5B zuH_x@F*ze2%P_1@d+3P1LVWArpV6nt$MhAdkf}YMMzx_&N3Omsu~$l@T>?873xeRw zU*4ONGRd59&`tI**8b6*64{*Op8uu~3cxHCJ6@V9>F|wn&ezaZ3zqa+%&zkeh4wd) zbqh;agMWX#`pQE`oVUIlBKv=?d>gL(Sbr?d_b$+mQ{{~dkveBm;;-%09gnpR;;u@0ncc0NLgiY7=kjcOO| zY;O*qtUJ|13(75~QkfheMw{m_tL8+kju0~OfByWgQ)B&D8U_)pv^YpstaF>B)WSap zcJfX=ZVpk|5)x^c6)lOIX!*U4$yj%belfO%JB(T+iD?#ff!bo;HiPB=1?U7`&J7hS zL#>6Lm`*w^xjCANf74w%K(%3rvvuB54W>4xZ^SHUDrG>M7_qv@XSLpR zSfoxzPgiU%bPhX>pr2{+((Z>hDt|N+bVolgA4EmIzcU zDHN4|sKyloBOEP5^d5PGOl;IGe+|&FYxD6&S8-Yl4W^|8GXoZ}ekwxzaL5CsQFZJV zkCnO8gAl;Jl{MQq`k6ugD1~!<7Y^YVgppp!Plz}t#}}c;-w>@Y%?h{m!?Drbsp^Q{ zZsbQ84vB_=d7R7kH(hjJ;+NJ5?d50MAZw&T-bS2Xl{kW8{P>udwc2IR$*j0#bid5D z`N(>nGuHo;^%GcPubz|s=DlN%`IN^?Y-)JUJpsV-O5Xxf98gJV!@ejWFTLVVu$FzV zj+6Mr3q&3&ymT*w4`Jsf4M(7A4y=7oR}Ng8&wOjG|_s1B&FyuFJO*=b`Zy zIs-bQLrKwCv$kLz_b`DrKplqBzIP^fSQ8VpM$TBG1VGgYnPsY!a=OjvC@OAQq(qgz zGhT4@rmnl~lkA0xCod61HV>&LQ8rc(D-ULq@0G50IZCh@PXwm}>uxtsIOemizxB-S zb_5J_dXG%ZDRRqqA7E1fO4_mJ{r2)j(4Rcp1PW}@-{uTM1p*GQ96P-qQLuxCk3t_P z)qX|vnxkH1S8au@QiT*AJm~I7kWCFEa)EJEHH^#gTn!1BzlQKjH zl?J8Ito?IsZh@T|V{&3Pfyz6hQ1(_ZsDXD-33?wevC);>A(U`I)gxtj_kOsQj#Qx_ zBuqzzE5LmA%55Ir;tMi=#3SEE-hO{A1V&y5dou^kdAZBm=&uL%1|K8Y=0X$I|n3Gp*{D+J7dU z+lbvpfw?PAq7S-`(@%=W_#v-uZ7uuip$*d+-R64+L_f^tBJ};%0%m3P%l%iyfD!p|y!R8?ZW$Az*pbt}zT_%yguDjdYiS)Xq%}*v zA{-?0!JZsTqqlA)ZnA3Q$1_@Ydtzd=NAQkcRGG#?Wr;^2qGu6ANC831b5ulu>epy(;k7UM)aOKrRmf zSQAjF+kCMeWCLzs3C{w$c(k=5j-hdcX{`Kjd?T&H=J}4H(hd$FSO;7reSMp~Y9RAM` z8tN`g;(0bMhvI=WWhAUy=CWHFcHC1XV_2Yy><6iw9=j0hK2MWpoz#a-4B+D(JJl0l4S{?OoQe+`(1g!b{X?RX0$ z15s!N^l@Ro2sRAyETJ~7u(7Wo!v3^ZWx2WtBl*!`FxI zRDft}a$er@D?oz62uw|36^SpT)<$2kNR12?+hLrUCRf5S5ggIZKcWunOC71V05XgZ zGX5Do)cz&9#>1?Q;MNlcOcrto{cKC^pArk%MCEdbyQr!uXjR-3;272Qs8oe}3Y2O`Cwov&v!h`xWIs)V}sbfGvUpk48>SKLVG9b+86VjprSkZEbQH zeR1eun!wY&vH+`(9BwsC(=FFTG-wsDzF`zNXQ9-U<7#I2YaRKiT9sOjBAeKSl36~J z2th+HT4qpy-miBs_3BgI#_J<`WQS%-ms9~U=2vp0UvF1PnBipjj&*4EdhX)_i!gGD zdWQD28>YEeN(;A*x*G%J9gn4sN7ThP=jM~ARm`f)H}**qRVN=Hym_%3+KQ#DTZp%;LV$5`{>3h*Kjla8T)?X9hgZO6!Jyvx795~b4-Dwv| z(*gU5iGA+Lj|=W$*L@D#9^2I&rj-z+^Vm3KB+WzHC74;O@>%d=#L8@MYH_g!qo|dA z0Z86X1#4N^ns5x>3|n^(TYL{H|Lev>A&0rbW~8VA5uVrA1{}aEyc2bGmiEchkg%mn-%bI#KhZ%kjxdqc3fb za3-_DD_2`hZr&`nYoZKqX@p{E$>?h%a%6%nybDRnJ_bFIz5`nK7b4pAb}P{wB!NH# zM^(z!93Qe(`}GvUwGP*iHF^$fq@Q=G@0BgTasNhj^xb?)4ML@BBWTxtin)=`j6}i7 zBYnLz`P0-5t1z?GJu&sB)8Y|O*?D$8y?L23WEW5x@^mjByZu!L+OqW)r-}*R2cu8F zrAV{%T>pNFoZm#M&~+X=;ADS^94@aK*4Zv9Wms|ORsIDb)i?5u+|~QFHr{H<+aO1A zu!!Y)LWF^)da`v)QL> zJW^JWOEQDJ8*IxVQZ=mg@&0lWVv^1Oc<{DkOg04-H<1qSp3)Z5_ zRG?c&UNRkj{oRwj7ZWcvK9rB(QC_JPivWp!iH6tO@)u099{L~t80rE_2+7|;Q5umk zPcl2Z(j?R>R>+=p^E4Eu?&OFOdysGM?HA259PkcnYq_42n>+mIs*T~^ zJMOo{yw*cM4ymKUOx;9w>k#-)t`hr3EY@&5mwaT<0Z{$J4fvaKlq1W$a z*80$w#tE4WNQRAi7c1iF*;+_GskaHd;lgF@FW)=p;%dkV-ea%r<7yRi1`^ zXnNRfc82HisQ;%IQL&nOD@HHH_qkK?qn(H_KkI#yP2m-HX);<$)j74OsE=t6gO5S5 zeO86m&f;5EteMQKbCO?|E_?J(+C7L-2fe%=WVg#1PhJt8vP0&|60{P0I5`G8#4S7@GDysOrwag=AblqjG`~ zG&VRa;lPy7^%J&Mn_gJoT#am#BCjJE|2dfCsX($dTX=nF=;T!`*hH8obX~^ zD@t36e8Wt(2JMqh+SPM5<0OV1tjN|2t5Sr4Yo`z5l%~) zd{@S-Xpk#>kw0FJ{bBkclnwS^N|(V+M6FoDlrg14exyf%@bd7;G1Y`5OlL=`^-2VS z@-9q?ks|}7H@*ppbYl-5J+ZPQ9r~g3`i*b+bn8Eian=s;v5^Wk(-xJ9+pVIe9haT< z>q||)*BEdq(%CVgK1l~`j0r8ceu{{Jq1m$+D>EX3;wno`zh*>wlm%UYm*ANVH5;{0|>$}cO6WUyo2!+g`<1{GN3v;Q}#pUFMt6B&xtO| zH%WtE1QHi`2eMD+FeG3MRQar$*DlZP?4q=2O9boXt7=pZm4I%yLMc^!#8I2Q%5X2A zJBZwm*jZT@|AFIYR9CtyOPV)UW)+V}9M4B@EmkkEwnpbt!#6N`JO7^c%ueMnn?|&b zX>elkdS$wq&&JDTpvfR+%qi5Jq}a+0Di!>_Uz4@t5$wSC(0*3|t8&f$l%9;xP6^iE z6jL>gk}c_gzc@D_os^jp`>bp^=tq^kS6zWV#UB%!1RfX-5D0RXcZkNPNPsqY4|F-Y zhKAChK$IAyQbuqC-nwbpxGyY zWy2eP%ZZfhOSkr>F9oy3{iPD*2M$*zLC4D$^om&2dqRj9oMe)qx$$m zicLq8FNw0baKHm_=22>*M6}CUk_btMuH}n0gMPohz3cdUj-f|awVwmtQ0&UFDWLjP zekDdX@!W^xBS?D8@!N3OvR62CJ2ph8)Uh}&ptaxQ0X4>Y*Lm2mj%`925tcHJ8B~g^ zFX&3i;hfVJ__y#Aog2)aboXXvBF*>AITgOku!g8QR&W#u`qG)*ljWux0IU~)KmOL> zVsy&s^{-z$_1bg&rifV<@ul2T*Of?E+=Gg%fA}I2@98I~^5;BlgJ8G2g8KdY@0Ek# zw2xCOYO{iRnMC(jg01Zo$kvAj>&&kh`-!4b@J4v z$szsUxwM1ZhC0~gS~K881xl2R%aUB^{5826mS>1ku(U3vAm)sH6t1+7C;jhT2f+bd z{$5rb`NqxRw9)nZ<-Ag@Rhel|V1IvO(p7Kb6Q5~8xBs5CtKVM2+9hA8(cNJAl`@SDLIe(WrTW4mOX`y`OL)o*c981Eb z8t`nNG9toK^DYn(!jvb;-{}B@(`t`nN}@0iZhVIdq$(lVPEgMH+kT0(98{2OEw3+u z9`~uuj@LsU*a<9{Jx5Sg_Yi=>-|Ivg4{`z)TWU?H4FpZ0M7AzCy=bM+Fd*uoPQ4Uv zBvbG+fC=Jsrr0oj?_>p`NH|2QNt-NH-ZUF7ZkI?d*(gY#VT=~s8HQs)q>`IeGu<~= zi{s*Bxqct}4MC^#U^J&x}fonaZk|1 z3}WSFqFOo2y=xPo$U~Ig^-2qF2A(1ceiJx>?9dAkP6;rnrCNA9aw(rVAibFO>gOMH zs}O1#0ENgamMNv3eg6fsesR+m`VGJTJ#`S}tIZy1@I&1F78m3x$)z=sN4xb5Y z5S|nc?=>sH_dczdL{4Wz+OJXWk?S6zJ1yJkY4UH*+s)ahoQz8EeHyyi?2eTDAD_Nq zrCs<1XEJTetHB{!``wyiWJYA~d9B-t#(o_;n1OQ317(2WZl;a5nDg!AO-<0nFMD4| zXFCY5>ohW?wgJp+g({@TLSa7|R~0FNu$@NfVy*fIXew`{zL(tZ=`8V`-dbq^rO?7^Z!oc03+*7iS;7|@SZjDAOI49ur zYhEkJvxSf4y`EwGEvRZL?-1T;#B8#cmtDS7)VE9B`rXv%_g5uwuEua2$JIT|v@fSw zMg_lf&_=o}em^wKn&eoZtISWP;sMC~qUzNYQ|-?r8L374a3)n$#hzGBDUm_BY5#bn zRXO}96X!hn4uq^441vjRk+6n_Wq`r%i0D<=tXMT9iX)- zC<6!pBsm`8hc@kYOGL$Wz9UK}e z!_9}P>U;W9f_=zf?*PNWFa zoXHZWZ<#c|+G-SCpVxQwMhtYJBmZ5SrfO2dFhZR($ExZlo{dwhfpOYa=g!*_RFL>w z9i%-m=r7QKY+K1)iu3`zm?5BIF>V-+1tT~$J(q#AGQfRjY`K}Id#q#2L6*9fM@~;i zGuFHXquT-xp)UYpdQciJ2eh7lR7Lz4^!P46H1GAaA7A1R=>|}CKt3oEIU`G9szxPa zU94IZwfHED$CU8V5xs2|bwb2Za^7v?z2Rhe&6|X2WWEwK<}cD%LvWa2_|A+H6K?pd zKP2{WBppde)y_JOt9#r961_5#F#vMiL)Qs1Fh*FLu>9L*RY|6AOc65q+Q0I(Vfa=3 zk*!HlpxZLog-aTOQl25mNT=F=E0#f~jvs05iBU&LnxJRu@aVwb(HcZgWa6%a`D2;z zA0W^K+q#xOND*M2+Rc8&lqYhEB@@d*f>4}!aCHhsGMlf+*_nV_dl9OpYdqVShewzJ zNc0?X=+#*O1-Ol{r8Q(S8DuAy8+#c~`wuL+N&z((x?nS_Aj>%Ktwd)+Dn>$)=h#K? z78CajAEN~uwVx8T4-1t_Md7<)l^+AI4+e9g@zf&m!$R5TCf7Y0gza8+9XU9}LJ>pF z%N=2y=`V3e&E0_$Y#^Pj3VqVw(MZgn8haoEoXzwl7qmitE5JZQhr>!XOOEJBa;wLp zSQ-%+h(Q{ELton&1W^XL20xFJB{(Kqxc}kgA13ZTci!qks7Qo}_7XIO0UQDCiqF5J z0IT~9?jaT&5#Tc_2moZ)q%xxc)G?uhp6wYKyeT_p_? zHFSUdA{^HcjgXSOYE`C9O-*g<3WuKiAS5!N*~W#LmggsxC?mQImF~IsKh`jMeHcl7 zJk$f3+Hii6c1&^JL$S7n`c6TUvlAR|XO z7J}j`wj4%Ru3hB6BW|j)13fC0$)>Og|9GHA@M`4Tp6$)-RaQjb0S-x=?qw45SXkk? zmDX0K?g-U+Ak4r-B)6X2(bmQeJw2^h?nEQ*(FwbI@ck%ca>v&NKuv<5L)$LvB`wGP z!7s(Xjbu&r8mLDta~@1b*BEEym1jW+gRl;#7D-0N#2TAhS@nQN%m=a_B-9M$U=6;mhPND> zrGtLdc=)i%8@G}#mzzyVh#}~^A==HNk~1ryNUj9>LKcm|Y@Ip1acnRIpDNcNZ2g>_ zxZ7lnR7I`;({cF|4@>R3eZ>Vu1u=K%!~w~K0_intFi;B|dQMK}u9q!)w(ut=fG+_t zd6mM~pj4&`Q|IL+Uu|7Z5ovg|xLS_Q6uhZLBu?G{-5a{uuu3rEEMTlrfM`Y0DV9}> zusweeGYd7!=ZFPG1^5Gomm0U`QfHMw(F|&74z(7e{u>W?Wf3<3znsNh3f2W~;SR;= zEi6H)>$P={v7BFF8+8tJCAN+guAo#1Z7?_{jIGc3SZ+QzFjnDj`eo=@bBqQMp;i8pEjOS&pA-wcHG^%@KNTrCuZotf`jYfblu^O znFR+nxqBBhM1B?=`^5Ij^8RNApS@g1&Av=Fq%<_Ryz>=r$#RNus{f^=bp80#V;q-o zl|Q2h!>*NQ!``%L4iNlyHl>CyttuW1a?jX<()CVIV0J=FGm+(7x^iSKEPKX-8}BwY zfDKknIvHos;@n!_uVEIXl++Hgq^N%#ZjRG!9*@fGW3Jo^cDJ3kn)QRfoH;f3q*Fq? z)8CK~+@P;$7GtX~lbZ%p8;eTuQN?c+_-~g6h7bZ-^PRhm4FAA1cHdC6bnd^Ecp;JJ~_9obhpp*`LxHH zo$tf*U_AGs%0e9r-?`GX))l|8Qc2b5@37nX{cx_YW2zMvUL>x+)@qA#sHJ9Svl2tA z=RU4|%{Au0Cn^4Me@6BRy)sL^t#qw|poqh;V|w5SnFyuS!1}@ytwDQY6hs5cwPtN! zK^avg_qkzqWM|GIGJ<$J{ZnnQv9J{$Gf)W_iN#&QWGP8T55pX@){reG)y;)D6z@=% z9c#8VH|xk-h#wFhiaCx#Wr^%}?fc9uA;Q98tocW5zG$Q{|E%Z&sagy*T1!_x(WEVT zB1)I%){i^HAJ5`WWfr?NB8nIK%J*aTEJ#W*?qu$4+iy{bk@M~cuYA==tr=1h%@`Cu zPsGh`JNM2@8Wd5A;eqwj*DwD_45UcE$NkB+YoFoSxxlwlQsIJ(n=W?5GZpJ$t}$^W z4Yq$%`}(MHykwX)JD(yC4{yO|N?pvTZ>=S|k@|~vep(dkN57OF#`_hH2%OzXlH~L0 zW$yc0i-9Kr5rjVbGlRpG)`ZTp`JdHfDfrpzCg*H(@^pHb^NlE;C^)o=G<#JhGEOB%TZ zGrHnmJD+mJKWI~qW_-F|-&JkiR~V44s1Irkt-%w=Gyw>tD<_Zv$;sU+CV8 z?UC&txg?fKm#c|qBpzj|eeG_Z2riE!;i$j3$2qL~&j=$xx96SW zMqA6|D{)q%GfSVRSs=7<3yV8M=e0iVt&JUL_KK}OdzIX@v;pxh*9BbTiq-_xdYHbn zj*lO6Yv0bR@sa)iJFgoVHIgLPR99} zY3JJ&;lP#leN*xzM7&L$5GA@5q;-U8TIZL#$oqGVo9LBt<0+Hnf*r2ghrp*@dx^H- z8j+Aq1Z6LkTy|tfRN?oZMD%=`FW!9IyX3tRI;AK`VGXg1O zM4XmoXDQox)q*>`++^$YOn+hi$(>F{$9l!`=dQ=Pd4pgnU;Rp=6A9{)FpTZytv7v3 z<=wp7D2S@^xo;^&6QkiBm^=SXm8+c9M~f10H!gSu@w1H029g07Y3EBGL^+v>W$n^q zi%o?kwk!U|Te&R5^cY9#?c)2uOwfOLcs0v3GTbtRrA0Y&PL>#j7E$qYr@hrE3U+h5 zl+o5cg>Ja)ovd2TL*7}m+_g<7AY=JW!shcJ9F-aGJfs+8xyfB_m>gI?n{hsXQt5yz zQ_8o7Z{GJW__wvyv3@b3LB2ubu1(3g2-&H(j50T_UB7SfE zTk6A#6*MG7A$5av*fZp7gv@z~{C^drb$y}N@J{=*ZDYkp2J0hGQC8jkl;VOT{Iu9w zrT4GxrD?9a7k4hI?42f2e1}+JCbN}WIOj8Oxbo*0vVt&&AMvBhP164X|B7@G{iEib zXFa)MS=AYgsL_P`^MxlAdG91eXkrt+(HhzX4;cBtVD2fYM)ONQ7m;2~FeJ5Q7eDLa zYF{{=ZS?INW~U^Iw4gr!?xPo*Ev$|qiPP)B*b7#&?UuIWjOknVGq&W4c7|DZ{X-~c z!Q|t;LMCS>K2DGuEu~xu>hlcVXe~x6FaxT-$i+st;+?p*M^o!3raOHERZ?q|IQU~` zF^dIMQla&*t}d8N-!r#QPn5;5R-DMD{-u~F=B;!Hk)#jVs z#=C|WdGtP$_kUImyNdBU7*2rubpKj)X4E@~1=q}io;i?w zt#B_74~ZjUfm9gj;qccL7c9f~_TR!}UKTn5@hAx7O}ZGEdag zvq;kWnXx%c1z^)pWmY_w+iHv{l3G}BVoj$Gk_!3*ny=$}h z0oYXBOPAoy))vyV-i?b6vMkOR!22>8uIzftws^v-HFLQCND29wl(RjbS4z;7J^6X= z2*!D(Wi15b;o{r31`fgER0oZPvk5_>eI_0)foyyY-NH|13wu#z@49~4VJjFrk4EJj zSI52-%Q6@9V)Q_Oz*pY}UEfub$8t&8{OA2ll7bs#(Q;8z1pz%1I4m_|ZCc$5XSMGF zn1C_3Vz1S?{sf38Zr9Afcuu-1u?^I(G?4(=lYuxZ8tSiN~R_N6Y0=4R|zzE7Yea zP@{L^&2XJQQDTCfnjXPFeTs|;HWn<(R2um5Ut_sGq*6VpQ@Evc61`m!dqhI8n8pXWzkaB=WT2AU)gfeH|Lt69 zJz5UgO!r7p)@UJGmdzjSYXh-i3>oT;4psJL^L6jnlC1mbT~&v`(nZ9xGYY& z!3fi;^dG%K^ULfqX1f7uZs;b*f;+YVr&ePW*o+^03j`N|f0Wj?cXI>xY3Cu_Q@(TE zvS(XZOjRag%51H-hRmI6i-Nb{2H37=y*)nLWg&S29J!gWj}`G2HE1xRoQ4duxm7-L zn9HFWb9HeX^p_ehEnV^%o^YUpEM6ZKPgnsbAR)b48)y4_$l#h8XIEhT?c8%aYks8@ z#IomF7EA>vz^b&*!6yBnbQWK25^_m#0{`n?+$nXUMKC6GAL3^0ce3r|vfCP1XWhHa zPli}CdK(W8ywKCfX$ueqY}GW zSZQ-v+K!lfCFcsa1O1(oqC_UYo%Qt1iAieQ$*V{v&2|@-X1Yn6BU((DrCavQycKir z{5ifs2ybPl$#TQaQ&ZZfr&F&Ow5~nW&R2$3BC4k`+b(_WRpxm=eklpMe2+WPDmvum zUMs(y#z9JbeyLZZIhbF_Ie4A%apg1{c5=xMbvfpaexJ(=Al%98-j0m4te-A!u`6|_ zDE1iEhkg)%`^{ImL2gTavo_4n?3XDm$mFMaCJKMcHBR%r(@+U7p2F&CXwa1_0n0J{ z>9u~f#xsWrQea~@w&U$qBl}x~b*7usRY!gwc=V5X3K1iy9ni=zzkLAZnIC);9l#V5 zfFt+HOboWR)ICf;mF{@5dc|zlL|k<%CF{IWc2gL>N04KPnGhsxN5+c3J$-c0HDd_x zTdwtX0|K&cpORz%*fHbq*LOd?`^0qv%$<%R0Z&slHkx=o%x`3~Binha>t2&xR6 zy@4L#Sq??#HaT9uh8!k6UWR0iLAK9$c22nQk;Oa-m$?C@0Nl3WCX@>tAsnz5z**IO zqQ4{MV0Q-VbZNZ5aLzw>=}bXEPHm zNC$a9obB42!8#8fqNrgSbLH}6Icd_oU=_&#KK!$q72Cw;I5FCp`_W@L*yDnezKl>t zKPS{sW!IAZ2<|e*hRG3>7Vkc3O-PRX`9@$ftEnDYaa+$$M&70mZD;y7b?){Vbc9|- z`bn0p_l*uJHKXVz2Sw$#*DIA3ewv;V4_`T6^1s?#;t^rnhKStjv_o{HyD$hcz9 zhnWQhVLxfI!It*>wyKL*ANFNvII-qtspt$^S_Rws@8I3pY=@!(L?!Pr9cLGLgjtL! zW^IWSYW+u0CH;1x*5{z2o!&ydGle=83ALt`Nxo? zhKludc4^G&Fs*^&epbPq@r)4O*`?a8GWZJx8htU28N^S;pucij*PJm z@(z<0zMzTI3;hyG0cfoQuE<>gp zhrPpUtyHiyv@3Y!#lDBa1=mt{Gb~wg^q{}uu)EauxTNwQvQtV&V9zq%apLLL560D+ z*AC#EJVWHcImJ@?G0|eh4-49+exIDE?k>V6F@cc1T<1Jm9?V^EQY-4)Qw>OZk$rwjGR_+bCKVzRM zCoCw4@@v^psgQHx*21F%CE^X(VQ)3E!F$Z=`vpZo63<}7OSGv;>-K$!DzMmSOk|{0 ztr9WpQaoKodFjiu0%> zG{=M@5lwe=UJCppDL1p9j_6if08ZsGOU&jFVI0QJW@4x|;myWwQ@2=VVzSRkxSZ5; znbY2lXUXdd$0*Tkc#nyvWt3*u4Zw*#S?+2+t{^VT6E!w(pO?sG@RL2`LY#RVr|E^o z_qf8*Yt9!d*Rl_7#u?$(8bv8Gs7MZ2$HkpB@-5PQwcSU2hAW5H_|Dd9AhI4q8SFM>wkrb${}m&a_No+LAFd%p8gvaz7FiXVkKg?AX0a-v^sZ0vRrakWcXv6e%c^nSZP zhf=viB}G)OM^J^-Zl$bIcGeqB#(e@JeR!MA)9OXT-xzrx2)d29?*9iJ@!d&dyA1L7% zDgMa*FBALD>`6~ZI(w|FUV*1Zo14F$)lYUP$2#r&zA9`9%qDKcy(=uAx8y%JQ#|1> zS;5^EsA-*WA@^bH%*>gSfWXyF_p5brX+YE*%Co+9u3R1P_&?K=Cvl|Fu(iGr|GGJ) zLW(a_`!Am&EWxIXBC~n~nFLGQOXtdW(h2Wk;sPOf?Wc@%2aBaH@G>?V3N$%G?GSq) z)iG{5RrDffn6C{HqusG=F0VA*lja3H*5+0x%>plnRZV9r&pXjOup<&-F< z$c7HsWIF~k{vEb?1>K_2TGQIv>zfS|gmFscbDox#{J2M{HHt%Jy!g#YG;9=sBI%W6 zMfX;_ChP4<>$=U6o$c+j=1%}nU~yY{>^3?2i^00bFipcRV7sY1+dU1B9Jsx8L2h1Z zyINGC23Jm$chOz0C~-Lnu-RFn%Wy{M`4#QsRU>P)0q#zW%b2Y{5u2>2QG;)N)z>A3 z2mNaunKgp%l!w=DzrVY40(C5-D}V~D8ND5g{ZtYxWSHW zP3z89YklVSOy(+rZ5-4(va;2?sZejSAWk=Y>FXz6{>zmUJ(i@S)_t^tyE_|KP!cS|+EVM-|bOg2HlM>1|Q~yFQhD%Q1D7^@!><>7=IPGgJt!SFcASkov!I?wsxER*4_?i?n1)UmHBv4kk1 zn~6p?p^VYGOc2?QGP%j)=b?gG^S(6RBGXdGLDW65(q$erwW>RUrS`1#s}sEe{lh^l zL7Cry?d|{>Bj8>?bFM(vc|2sk2xRih!Bf94T_{P#oyj(7jb58G!7zfMQMutmueje} zw`XF7NGwF@lP<(~-Spc@S~Y{#;l7kI4<`$z4)`QlXOr$%(=58vcblHtoz1XG$ll>la6RqdAY~eQESZ|0;EoA32VOxppe|$-)s2*R$tjGe;CRgy_|1uoWtIrhZ@| zDk#~8^&2K5^On}GTI=DOzERdf!IRNDl=w{o;qtCyhzwaMC;3@x7dk$LA|C;B#)S?k z);Z@g8#bG{Et_Qc&l`FOfTiEpWZGUbR3=I4iJdz-Z`JR?tvYi*HX88$@@(^=Tp<*H zGrz~&_w8r>EZa-X)RQP1L1g7BHxpP%CL_{J%~-bZeAz7wZgrqRb zyc03zP8rdi=%eU1PD{-tI&MEAv9he+Teo(ppn+4#{*<3yUavYNy{Ho#O;-9wK)l7) zR_l)2PQw69>p-04K`-Fu-`kP2aPF5U2VJ^SPo^szS@Esk#FvRWT9Jm1RmpKc0-t(| zo9|#plYuB?Osv9Vvh#&XMM-LaWO#nJ1>%8+B5M5>%SPK#FRYS7Od4gUSNDp2V7#{& z(TyNEp+`-{slspmAVJtCDHg6mqTubhh=fS1ZhgKb7us8+5Hy!{)Y9v9j$~fpt&Ky z08h7YtCO^;C(aCch*{7(M+9T@uXEKc)I#Rg^soQ3D)6C&*WSP?nUS-xv)QM}V(?*? zIfTh$UMeR!4p{T0#RsOeZ_Q+Go9)cGCKOk$npdj%E1Rd#f2wdg(s&ux*|qWd1MGu6 zG0T*AGVIXoG7PX3)sAC5RAW?XIr8_*;U;eqc817041M=;>KnP7d3LeO{ZAi}y^(Y@;fvCd+hAI=uoEk`f7@>@K(^7p?g-AK-(xlr zuzy9Iy`FZ3<+t6cVg*`B*Co({eAQ^qzQ7fKb4TBrhuwp8o+}|QuW|nTlrj?;!Zi_} zBFcX`KU-qQ1i$Dd`Wd!_*?R}+txrG1TWLBj!{#4ReZt{BLhP00&ac9HTqme|{84Y% z@j9A7wnu|O({EV2cnoH$Bygxw#GMyMomBsk#hDzDKx2jiUtWjAb>n`sBa4ItedSKK zp5}5Csi`{&lbPmf~kW;M_skxuJxb{(&P{Q+ur46|(?z zdCUmaS`=$82{*KE0oZyF@1o?eR0dw~JZ5KO*!GV<%T_TY%$_S3VPV4o-*HQ2goyO$ zX|~f3g^Q>jy(2a8M}!FpRkBA$;*Us9Z;%7J6u1w(;5*h4kFk3)n0ho}onCD()`Kzq zXadph30pMcYB3wSocHXXL}azPO$w11D@lxh#K2YHtT%L%OIqJtj10M?I_(mR30z26 zeL_fB?0ED23RtOS{KM8)t%(m?KZ=PA`)ekgO1EM@R6tXY*~~6SAn8Bv!Dk->o$-Sm z+T8hLvUFH)2CJvpmU#t0>+ZRK-8lbu0|8 zT|7qr96I3HA~!NyrgDJr<6#@weG$@!49^EC$3EXgPginNU6sDUYj-JjB;tozsQRKq zHn(tjb8ohz;Jr5h;LV}PS3~reeg-`-@4LJ{3c58$w6{pnqI2G}P60fMNZW4|cAJgA z0&lHQQpi`e``y?}tG=QJ-xku3RBgy!3|F~9FbL4wNYI>07Jtw93^=PbyC+{x%S!0wbL4kmv;S?cMUrP+9#1|MVp7lDBP&u^H|br}%{hgx zS(8m-_>!Y<_jpXDk#*UOH*X^?Y%7L)8tt7dG@qCECH_VeT@=?F5)dSxG^i(b3dtJ4 za7A^>g}3{)R_pp6=41AuM3U6YXH<=}GB2yP-{;NK^@m(f57+LgkaFbncJ%4F`yE@a z4$%Zpgj{;_JfmbLWm}`>EoN(MH?zcSZY!l`s+3Qag)Et-oP!-U{~i1U=gGJ4h1yTK6G#)4?QbCRk>ItHKcoP+-77-G_Lhz#d3mFnUd+~S zju>;FAg=vUG}~q8`Ptl-VbX~r(}>VF^jFTu=pt25w)5!JO(;$5m)?I@;%yrq~ZywQT0=ZRXP-ABfLtfsB&$vfEKf z7;zxn$bMa0tH8&vuO6zu_yV{X@H_UFy{=gcmJmf-%B4uuOSI+2Bt)@4UpfI>9usp4Z+i<-p%qjob zpWHD7=W4}KbLm6&-U3x5YlxfRcBwsfWAV$@_qt!gvIjYtB{Jchd)wD$XgiI*B?S`N z(^SXNnS)G#5%_<5hGG`@H` zF`_|Mg3~7Or6|k#PnmJmOsr}dIf}?rjj?te9IQxPYH{K0vXh1|z7jvNxt;(HX8W1O zhF#2B4abHF|H3|&wJJ*bP+&#C2GYkO??UvYFzMa|yxvgn(baeh~(+2{^4I%n5Cf;$-DaYPoDBSu+RlW2AMwT$ZkO z(UK@uBPWS&e^$3)Ocy1>1!3VX@IQWqIToj0QJp<&t=9PbJIW^c{*BNsk=ekQt*4M> z5#J43!|;lUM#@y+&8v#i>rXWfaJ+fJ_xXp|v}Ve0yRna&SafB2Z*RqF{?aBIXCZlR@x+fWWu`~gAq*FmK>$sQ{Q7Q$>N*apJ3QmnEZa8y=PVxb0y_E%njo#`rOr+?V||5FN$oMk5(NR> z2x6meD4SY>TaOYVSyy7E&bYRg2g%;<+1w~Cvt`PEKyopAorR#0=y@D6{cd9PeB;>2oFTgIcoeZuNoajZsd@$k5ke9926t6isF$JuviJEn%m90*i*H7pa z5_cth`fj`ZFw}42a~NN<=T6Zut5vP}eUJC<4E`&88fBx)Rk=O&zIbiII@-5iP4wKK z4ULnTWfnXh+xex2MUQw1aFHzWdk+Kz&J4KNFJ0So}Hnm2ZRii zxn=|4F}g*Bd*AOYH`Chgw6A|t&%N4}3z%``aw24FtH!u9a4n?SxbNc1Ud+{D*;;Ia z$Ef80WA8n~qDr>5VP-t)jEO->f}nzkWRWOwRH6t-YH|<+Y@&c<5Y!O~0wMwm(u_pO zvB^m#2na~dIW!r`AcF6z-E*FE-tWivKJS$Arcf36glJijQ2I1;aG?j8vajYO zNdpe=+YFJTM>)MM6P>rjIdyS`6a8Z9!7IstUmOc;^# z&u=a>^)+k!J|=ey)!RK_2E^+XzWv%Nf8-2}Y`M+rkMyH}c19xgTnUQp%Dcv0>S^&v z32QX^*~wK@D{h=yaqK%hHpksuR@F!1w$x68YGn4Zd#h)owr~*&QUBb0!*2Q>7XWBA z%)wSoYo&`Npvl#gmjCY<$rCP~L8)&GHb16g<;T|t9f|Y6`#M{8$NEP*JY_p4*n#3$ z$p2F^5lA8`M%3iJb%AI_3!A($^o@eo5XK4F4g1<5wGDkPK((>41Rhn`wtRLU8v;e{ z2fQGd3^|v|g1sx>xRJ;tUsgY#Yaic+A|nuUBM3lbZp-y-Xo?D-Y{g4ialY|J+ZzDg z4%3v?XQHHc1I~YW>bQc7E8$>uo_HwK&Zs6NS8id@)ByyJ!dI!|I#=k=_2BT8%uo4h z>2lo8YrR)-p+1|ys9W1_!}JmF!N-yrWuP;kweoJXvM>-cBv9j{>KfEkKUZSvfm0IA zIbc6kv4%FyT`lo(9Og7Ag_uyz^Ag9Sjb@~J0b`tWGk<%xX!&{3&}Ci3vlwlMut=%6mrL%kt;DB6@AQyavCo8BEf!29=9Up&2TJcx zH{6?&h26S@YALwKCfK#oRKGKg(*W)_2!xl`#(#i?yj!vYbQJyQu2q0vV8{=qUTHlb zZ!w|a8uvW)som+hSz6ih{WPAHquC=&!O`XhXiy?l|H?=Qb9m(qj=CbdKnLt{+TDwO zWIbZCChiC;c{_%E{9C7bpGSHOorB5lZ4_z<&P>xs&9|fGGgX*F<8&^iZ%kI0Ddp2b z+zAD?K_p9J3AVNDRW1^hdkbEI))Dty!H0FcwBHxH5XxlM-GI6Yfy0b+&0JE=gY+|` zij|q$#0+%o6)f!`2uezqEk79Cn7nm*PSM)@6wta_=Ur-@J^iqiRQw7;VOg-rJP4YS z8V{4b71AAVv{4AwJcWsqwR7OM>egF-X-PLRzwhq|Fk?Rk;pTKoo7q;7h*5R+YadEI z#wV}$2}lsf8EiEp=wMHvFzHaTzGLh*Y4fAUJ}%&$2ixj{cO<6lL1z0hhby2B2+%p8 z_RL@8(K#dLeJ#iGJ%c!en&8u0qy86ZR6t;5xyog{02l2NfV(!eXsxy!j)Edf1m)K1W0m;-X4bq zL;s9JO}7R?^l7ew_n8;OoN_(O|5+q4O(i=nUUT#H6A(>P%f8z-H|a1z2-_~{SqjVq zu|iw!Hkfi5F-dk*0zvL28mX58>C?30b+;n1L0eNaERyp626Z^@N-+>2x?&!FV+0lPtQbtLM(%knwToSn88)C)Q|<*3xMwM^&XBr9}pj8cE+X79?*cJ@|mKnYGH zZ+PdnWKi70sgUh9>!L^}_{r#Ds)aj4Z(rkqGZiTuIdETjj`ak5Uwtt@r!%;9(DMr6 zjb5OKP^|}xU9Z4N`zi;PJewV~5+45Xfw{YROhBn^i8hRm3KmGQMbCWzZsB9&m3pwp zr*iAA*WOnfFhwDr=B>}T;*>qar0-|z2)#yuk{mmgF9iwV>XMQ`BhdHbwm`S$Tm{EB zmNG85%5Bs}lT)p5v>*1dIL|#h#zDC3)&MJm4wxu9bJSCh<&w7a(|QG@3=9#00a#Kx z{f(QGsXN;X-YI4-sf@u^N#&bnBDZQFw5i?#%FH+FIt*4=|9vnYhA#!d9#U2}4-heY zOrFFL?Ca0^tlgFHGCd!?aq~-2*-9E$QH|Ulzo4wCj{75N^S48Oe-Vm)3fP=AS~_)m=k6QluLkImuTcK>MTr# znNtyO(w5aTs%7e-xjvfR(OtbW0Aa_d#s>XdUF--uvi87}PL83=tv+*L!momreA6e8%4Q}&6`talIr=oR zapTdkpAd5sz3-V=y+G;GT^ZmouV>4YPkj&A&K*W?N~e_bOsJ=M?Q;h%nE)ohYZrvn z0>7|hiF)O1r8^7ZYEPgZ@#LZ%RnY=tfX}x%rnU@y4{$t2$9EnyM0?Jyn}YS0@X{Ps zIsa8oFddN5|8g_advjmIp;bd53|9ou0w5UU%fp2N7v<)f!r>>QEZ93Dqire>2{Wtn z2Z&aHAs$kY=e#V?vC_BOQe*1r9-FS)I&MFb$75IJl_7ytR|L>7qQwPzz@sn<$>$Dg zDiQ$>DST}HC7(S(_p~ z!9ELrRjS z)N5xUXw~=vYx!7wygvx;o9&>(De>LeUjz)Y>!dw~T7bK)iFTuED-c5ESl~3Ov-@P3 zw!)V3kEWnt0&gEV9!f)53OhmvF1r9bo>x;xOy?Z=cEUjxSE}=aMMl@If zJCht3aIfiMjg@uadRl-IXBmFb?1)uE*HHqs@TM=xj=NxOO2m*BNI+Ct&7YAJVOY%p>A;YXGec-;Iu0Ik2x=>JG?H0~jfj9#fxA z`Mg#?JMI#fDq0&J;Ta~%Uy4?BW|&evr?UqmD^xg{BSgDn7k;^B@ zitw4%-{aZew-qRp)TO{h&^v%SB;VvS^WDCdooUN$%_WA>_*)XTJ4D%_A0#=nM@%_osRMDcZ9Qw7PC$z@vAr8&BX@v@~ zSC68*cnZt-%Q|^nWoTrnz_+_HgQl4xox@ge7mUA&H~C6x*%niN`PQzXH6xYg`aVdx z9uuBDf*bmxE`qV|p-g^j_n4kw8r}L7wwusIcE4j56JtcFMX_EbkJm%m8(RrHqL{&LdD*bLuZh zA1|vO`JQ;?GwNUtLwkP*-Px;*lxGN`9&B8#*~3Bxg`;^`1Hm#m&$OZ*ODqnDO;$wh z4fw?;(;%fnkv_S-uKuIewWaROQn{wFbHjXtufcsHJ#%&IrtBbWER8w| z{Z-g_aB*poc2>$anJBL~gkB_y4ie3TRSVcFNU2FNxnauk86qn&Ave(h=8ExuXLs+Ife-kM@yZ5DwOZlHDIS= z=ll-Bon{eL8n#iX&X_gk!RChu)nMTqJVVfXE?u?!;6Y%#foC**R-Jy<0E0fT@{w-K zSy9x@NrUAVwlQ!&SSX!;@LMF3)X{NJF+PspG~dn^`Kr#V78&QVJzVb(!1PP+nM7_K zp`b+xmF9w`+6br;Af_2$F{F=8=vkvSac0pgj9Zf8bh1V=MXn2H((5tnG4Uw`QYm_6 z(9L)?^VIMyIS7+cA+%UH4JA@Rv(N6^RWij%xm{^#agEn&Fnw`h|Fd-vW& zjRJtv5wxgq@P$S#At21^`H58WrcOzEBXATmJUH8QX%XjF+B&RQW`)N# z@3u9tUX)*ai!|Hl6qiSR5#ddM{h8=s51y`Pyr%U!k`KQ}hM)Q2RPwN#rP{(Ly!xK#xEuumfCd9t#q4D z5rvw;l-fWDGlX~ux>e1p#MwTrd#%&vxgcj7qFI<;8h+5_rfGWtX51tkX>p07oOaKK zHI~(aBbdJ8`vm*t(nCivi%2oZ?3!^2+MXRKtgmc%M;3%~&^0^yCd1*!POUmOqzBt>H|SVm&h(Onzk{F7m|d`O=vR! z_U3vBe?z)j4S4dIDpa6D0_I&Pw=B_})96mcAU@Jmg#!c~VG5E&B?kh3LI&~-t!#+{ zx3jxXle4~AYxg^BJ>CBu((<$tYa@?70SgevTZ!(0RKoe2lxJdt>Q%YrfF=zLx$G({ zrDb}u2dMPq{lRp|Vg|~vI8llW%z9uhc6HDjqkHwnQf?R9$JGHK8P2Gir+i)?X_1Lp zL%dMTI3tjH8XShmiHxnT+1M0p^oAdq@=WfzOe+pf`%_q4ao7vjH939i$j2){{W(?tRk)aGyVVgFBCV zTpqObzqjtucHkHVf)xjcB`Nd~;xCzIs0C*>e_vHb_s%F<%TqGN1ldaoF-APxAR#Dt zM(C89B*LpGcpYF)iaR|Ziz-D#k9R=hkXIw-9<0+)ZAl27dZxl72gxON%ntR+;SPP9 z@AeT^3VByFQJWe`c?~3*JwDWmLi^3OC-RmP)~a7*@(Gcbz)w*RDl(EjErh+;C=`dt z2wa=zOI7_2X*!>XbE}vg!|Og|N98Jg^^~eAPo%ZHQpMWYiia$Pp?DJ7HRvG8zt(w; zMWwZ{$pM%aNxH$1gMjG-UnPA|THkU@w!eo9v~k~%1ZSI{+P>6EVor-%+XY2dk|0FAX@zd&NtkX zhAZ=9ars<5b zYU?5XU}6-8XzpS#N-?!EPr(EqvxL6guhy=s#ij$d%Xi3W`w)MF&`+Sa*ic~k$isBK zhIlpLnW47TBZPgt$Hr1$8#XqrL3&5zUV_6OAq<=J?GI>b;LX*q^7K=;WoVl!cNJ}EI80i_|w1HF+Dh@0vUz}}=<361DY@!aXiRma?mP{HM92_xnME6MB8eV+pg z`HZHxRRlL@1Z=#)HMUXzyd)fHCoZWGGuoj=u3&^kOEtEn?G5Yf5rkD55iqpzWKy99$YT&9t zmF0!3;j85+5sMrbP)o5CGtXhq9EBWOsQy8Q+={7DhW)qMzSl+zGQ^X}L17ZPC8LZ= z+)}vp!Mv;C$CAHntk>b&EF3k-jj4yFZJhD)xX669@(&0iNK|d)wH2F>+{3)cZ;ogZ zpy&tn;;V(((`hLtsuhr34aV&bu~c?&ZQTagHBmGkoVv93kRI>H*iVG3)s(P((g@{u z21xKV4(X)Q^tg4zBXF-e>F*4%k{%4Mb>0%3Pt{YU0)!^e+T}QANl@^?X2o@}YxYg_ zyo_s!i7MnRy{7LDMB_iIp7XVFf>BkuzVkZ!Tm<*31Qdqo&x%J|fx*9E995|JM_QVo zF-GkQJv9p?2BRlmK#?b(yBC4UfDE&#{;$Eia^uMufd1nf3M!k%Evjz~#p7vfv97#!kL!MgQ&*Gmp$3y@| zo<^{DvHQpjW{{7@! zcq4Ev=Ok(lLV6GQW-7I2i!hkL#3|YcG{RYFx~`it?0lxJ-EsyTDArEF1q270@Qo)D zq&115JU+~Ik5TLMo4>`y?B-HwKnE#(#SY8&20M|RfFbySc-Pg>LGcEiskl6Vrpd9S zoY?|#*}*+HBgsj-e}o7gq{!mwfOn;A@^G6Dhm?7~sY**>!jNHBd&hh8Tr;JAiyz2c zZPZzDZJfOkgkpvWV^FSRg}jCX!M?2JCW0q2g)_iM;Zn|GuC$4hqP4WKkE2**f785l zs>hm^xSyG$=0Qu;k|Jz$uRE10j%}=m+_eAt-W@B)WFr@}Xa|A4TMq(R#N4So>lhsd z%+0W9YxT{ra+BJ`mrd?fDIP~Q7TPf}HR(`{DP@+b`~%Vk8lC?xXqzt``m_QkYoj)h zmGByN7V!Ta`Pn0knVHgfN8kd5$=W3rq|Fd83Kcx z%P@}{*ncZ?T=ns~MPGh{%jqB3?Ql3M@I89lJ29ZoBp>fwy7GDDsg4Z^e;U|8S`$QV zRQMi#40#q2;W#G`(0vSA_B`nE%@?%hfl-nwD%{1{!8t9SUNqL3pnV(zspn?B&A})_lu;P-w1PaRnRTzhy_#Dnp>S#1 zVPd!qqS*!5j;3=8y>P)yxH;9vKgUNi;!xze7S)H!iAl`bHgFOsT(xgCom*yyt2 zT?31kzuA<0)}WRJPD-)W?tM+75wzzk1io;fLPd8g-(sQ=vpqL-!8aIhRY~J2>A;&7 zafwUNxFX@S_69P{O%@xBb66UW>+juFwW{rEh^mJ#DQwQ*4M{5+K@O&+8g`RG_F;8HHtr;dCN$@RF35Mq-FQSp%w`H<$Qr*AqgX3cLtuSw57Lk8b7pxy*drY~^5 zoaH(@FZ|dyn2nF+gp|-JC#5DXafprYvNQ4z=KO^e$$V{!BTgnEQg*|uIAu#uc@IMX zdoSz7+IdZ_IJE_{l)6*x77C9S?){@XJwQXto*`J(LPlhhU=Nq>LSGsN|5nR|F=R+T zv^o#efu=R;g+2j(YLP`IR6>fJtPuEGe?boT8@_o8wI0zzWJPp5!GF6wfn@qFJZd7> zg2!%1oPkKXDTwse{I&0oX~H&!70tIYht;U5g78$R%f83Y6URIzd#@LW1NDc+C8TRy zNZSs}V#Z=rTr99^G~!cSfgAu5OE$0Q@8T6nF{)6tYzL$=!-iPtns<8l;xi0Dp_oTI zB2N1}BAEs%)0bsId6wMi3$3Fsdz+KDVlUttu;v~!nIPp^r|CR{=iSN;9e z!dSdAp9|7n5jSax!|%uj&8U)W)LH2gdiE4?wPu!fsZ1_cVaS*N)N~5&@^uDDm(`z} zY-@Q#}t~3V_>N9>ajIda(&hy zZliB5e@mX$U)f%HdXMDpE@hAnWmw3gbD92Mp2hdTuie{-15fPKC*u3h7zgB7M!dfp zR6&st?+TP7zNB>?|6akG>l05VfE@c3TIdIQ;1^=)R1C*KkHX;2Mqf&32@XaFo$p|p z%+%HC$JI2>1N!Y<>C|<&LVrOj)v1kE4@KQe^0^v3Z=l~F>Va6I=mU@3V;eHm{QP=B z)~>pdMMcnvrkoWA3hvGhH)T{qggwvvkgc~W5P^YbNZQdogCI>UF zF5makBgOL+lmUUrJQ9#~27BC;2i8h&zD6d%Wg%9PvRKK31pX_UiZP?hRm!O*k&iU2K6mpP z0g7_TinKP>zO`B|HjURU#G%dL5Fof*EHrb9N`DU2ThTtj9{)A91SfI-E-wT%solXH zxmQQ{-4-b@lNZSLQYWxSIM8BnKRq{IqkBPJBn0vaV)AIk+#yo?KXr#F4!wkjvp2zS9i&fwQkSj130zwbEizlFq|@6o+b|_C;CT z>yQlfSJonBbqo$q$zH*En{0)3v|7i8uW0OMSuX6#oA%Z<4HpGGR z;0tvNt!B?htepic;2Bw0;4S2Ml(ubp|uk zpI9QV7-ZO7T>O4BNEQ`a1BK%^A%IpFkEdf2F=kPp!J$tWdx}W{>MhVbqJwX0zQTr8 z$Ccyr8W=VAHbhJzTI0|`2z#$rj=F%dXlTN_|LTv8HxNXEZQ?y8OJ(eBDh9^E)nUNs zAO(<&k}yJ8mn&Jlp#aK+A)URJX?LsHxv^am&)N-zu(JBtu&@XtaL!TxY_4#8c7j?3 zP)`NIcgMQ89J*S@-s>EPRqK1aR^AB&ncqy!k9GVQJox4*C3AX}_{Kx&Z1 zI|l);L5@wydw|u?O80IB+FK~*Q~kkO#OR8$g3ma?KOEafULJbyvG^~1&Zeoc^1!ly zFwnb*s9(;ad%ETIqTOd^yr7;9klITuocbnj10Obu9bIHZI**O#SkmwU{=+4$xzrxc zci1~ri*d3nNeME6O#-j*ypzkJX}E14fu<&`RM1x9oJ8$G8V|9da!BJsLr3l4RZ0`I zIOQB|k$JLci4G0EZ?FQW3P+DZCV`+Kfa5tjoN716-ls6MqP8=F;9aJomD94`;_xod zd(4(=zGASf$qm#3qGni;T6w60g*nXwEYTg^Mrj#dBT)B`zavDCjZ4I2D9^sYCIu;J z8({evyebgYRHpGzSIQc8-+D~u<78CStMb%zUQfx$dHYdPdzqeB~l8*$JDqc(Irm16DQ#xCWQw_hIrI3U*{K_M&9#p`rN*+L&2<6LysEgGP%aUCb(8zcLAHtz?n&?+w{@^i@fJ(4LDdx_&!f3AFJ$y-#5tHwL?a4%8HKoVLBu^7S*6Sui5)#s#EI^-{5Y*!!Q;-|HH zR37V9W=H8hhSAsQQgbmd#_Q-EPUdy+&Jm_l;LtjLD4 z8theJ>KQ3e-6qI`82^u~qEcohPgOMYs24cK>m(ES`>_*DeS@KDNKnBUsH`Ci{)u_@ zVbn*Wz{ngFKATQZZ;qewRzB77_?*7h`>o;vXkwl<3M%`S{KPIUPU=1;KL;Pi$bu&vV~hV^QA01-*}+yja+9P z%+~8YhOn5~=0qwCS?X7;4cBgO(A>PtoB^ZmG_VReoKP)9%BxI#EsLj8r~7Y(aU+j+t8nk8GRo!BHZx zQmQS+G^Eh}nySr3Cgi2tgHtjAGm1k4xY(->lzp=5vOqBF?1QyS)MVXzDnDq^ZL-%( zGmpHb9j~qY9Rp2zZ-F2xG_ipKiRWH3P~<}VN{WoMf}*(^8*xIc71VY!P1UGDSWZtd zHp?~_k)=#c;!rPPX2Z4&oMD(eykxFHV1s8BXbJwS>f|_p@KG=y+8tk{=rza3(zWfP z*g!RBjJ@;tzNKZ;5DqGH?wz_}j#2MSoKUQxlQZJ;gK?zEbey4T0@z@j>6cPyXTc2^ zuil&{tECR)-X!L=#!B=AAf(SEB55khtYcAB97?Oq^nrhR!&%82F68uW7qrU00sFZj zFK^N*9~W#T+kg#om*yl@Xx!`2cpbZ-O4~%{-J%m@qY9-fNf(qP*?r-5VKH~Nsu)rK zea22t8B?xZ1!6G$4+*N5e2=CNywg!Wb@(P}S|*F!Vt=&ob$z^=ruZP&MeI7>-F7H5DaOY$&j%AgpGcPaCBw zyo)#qzRea|8uXbsuBEc&%2nsg2W4Q_j6E~+@F?|p+I6Ljrd6k_(C|9T-k#UDA!=0E(sKbOmJw0aT1a$P{dE9=C5yg&WD503b#v+|z z0WTcmFko90F=6jLuWufc(e;XRBe80tkt_Q-leo)zcAlOEtr&m)Ydr^^t?_*j)zY-Bq~K!3wT z^+>;g4$;!VA@@vnB2+){UY_jMx98zHnEPGbkM!}3>X99T-e?T6jtIjU5G087`9bcQ z9nrUzy>=r#7W%|AF_Tg_MT>ItLX@G4uDbz!wXS07$*+C+mSztg2-+2o`R!}r^hJ!N zUKSKXn8y%XQodccWh5NEq^2WFr&Q9%s1EDhRTQjT2(QtO? z<(sr_Z@U?IE)u7<8ArC|c1lk|y8lzp7VU$kiv^{!%b%NS6BAblF=!w{liK}EoSkf$ zv6Oq$eW9#=euEH9Sgo5MODd-Ici0F|bf#*oKoP5r?EP=*Q*m%&OL0bfc%_In7tP$& z5%=^Xy?27b5 zd~^;|GIJYiu>Gw`<-eh}UJ1E5YT0R+$(S=9=#MhJXZutj_3hXelv4My`vm!dkhYh( zQ%Fe9t`xh-&fc`C3+Gw+{G4lYh1DSX@Q0i3FA7A0v*!A+j4rLdJQqm&@BN8!_;zKou%45X)718qLd3(${f+FX2}mA}@YRvVyY{M= z`9BPDg{CQRUroyxNz+j`BpV{fte{@Uy&dFn_I7r&^@ek9CN@fP{2i|6SYbdZwQZui zLi(|Fhsm4@syXWk1*H`(p+|6~%`mAFAc%aw3n~XajCjW!feL1iw&&QXJztA(QCG@L zkPqF=t}LdWK}RoYgtTh!-(lL&cYsvC^@rFnzYO$GAVr#5z0=70kT~Ht#Hv>G327eQ z-{%W|{c0oJbYFE96y$sgLdg*&kuT$_UPx}6z_t4)zn3?*;J7re?c zy2B%lvwFL^^%Ics9WzoY$DdWP*-wA=Y-qk^Yn@rBWOvo~C$XXq@fqJb?g@^=a#=Js zuSk+j0*6du8QY9aWT|CqI=zZ)8X}HSYK`aF-oSVZCi(h3P-?%{xzef?jvg~E?DIGe z2TZE&L|>Oo9j82SjyQg^fiQI*+NI!{`P^yx_AOpEvufWF>a{bFF=aNhA?;sy%NFuh z?-BLPi%EywL1{KMN~P}M($-OjDPnYi-2yArYE_kAi_;%kY1T_>Z|$cM7BTI z82!p_~&6C2-Jg$V}waqOB}}Gs8&U`$j}3g=OH1H zW`~qG_YH6hY6ltMDJHUsg(Tn$K-KHALCuIQUV{p6Q8@;TM}D8<^!sT=ZDVieTPiA% zwnZ)CXaO~lreo}vjPH!zrOLFUDuEPz_;4loi+W_yv=R*G4toX`n|cb!`C*F--t6`V z!sKy#w;Mj|B1T)KJrOkC^ndC?4l`h+8~gU1Rmk|NmKB#83UT4=>d=u#Z0UALl;cwAU!M3~2g0z^m&(jwdfO9| z>d3mu3CxR6k!HR$WLG4?5G?oe$P)|5Bc?kmNe5Z8wV?*$sLbvxfbI0hQf(1JKM@LN z1Hgv(FqA0`xJ<|20CfC3gH$Nv3)-+);i z*!T2*l2Z`SAD5x{1@hlKa0y^zWhg zUnxQ+H1;`m?*HpT+x}VPx;(0`Is2<$jN%Eo!vfDkwKP0oA_+cmBry@)W+Xb>E}sBhd#pR9cpn?&-6<#yvA^7^PL@ zoXx5D@&f+xnO;`-|Nr?<`NaXx|NN!+GVymulhFVELh#y&g|Pqr;^OB=lK=Gp#rA() z&5{4_FQU&K-c?p#-Rp=E$Z&nTzfImh*i~-V<%tp}bNbxI{YMG0iHmVnUSas~xJ+%~E^o z9XsL9p=YjMsT!M}X-nwJetCj(jZ#rVrc~zpzZ!X>_>2iYk$<{J#x!=suw3e47I(4P z*L?}=Y_>Z7kG)zt<1f-x8{PW(N>RMc(yrGM#^8AFbbC*`;qt3x_s6bp*czY}mTEZc z>oQ&w(Dz2baypJpv4bWq@@C4y-tJbqeY?$m^g^ige)0Fh;TXFM%FlM5C356387o(Q zpG)tHv9Hwi7Fd2+FtVKXlWq3TJ2h(&T-wtI*?v9#$osOoSXI-rkPprkeo8Md5S84e zRR>wW>B8}9gG5fDl`~<<6EL{>Jw~3k!w!4xN_?5lFJkd%-UM;Rs|Yj zdvndY*zpEUgw}zNd*>8i&PjlD?*G+HUHN0RX)el6YJ2oaUx^dRB1I)O`WZ+Lxg&bN z9NxT_RD58!F@h_9qCT|s1d~_+Upd^fc4)G&hwr}S=eVCgNTlPWUnn_aVi)POmgi{Oq^RZr+F9@Ht0aU2DJJvN_e%((|Eow=$k(gy6=0 zoW^l;7gj`pyuK z+6|ugbvJk|VQDbp0{+%yMY#%>(1Oo;n>@U^LGk@N#-H}5nxZV(RS6jnBkG6?JLNos zlrZ9#vO!$DF1Lb6k;cBf5$NK>rYVm#^D~?%?Q-I*RQV@-7cH-GW%*hYUzMiQ($@NH z2VG>o{_hzNGqt1bi+asCCR+DQ`>;xBXB!P|>2u_Sz*&W@q}e%AoDYGoY<~az zTPbf3QbRa<*U9LqHACppzDRci=|_>esQc=q*T@DjZGvfPEz?8a$a|w}+{<=sx;dun zZTe)(!GKvJ(CM0Q&B?3ki} zVch=5Q%Y`)inFsKP|h?K4tmb~UigfCr9kq-g)E<<@vBt=dTrri>AJfedc@Xr0K9z6TD)U8XpH=Z~7MbZtwhtrWZWW34qm zMK!MFY(?XUq9lhlViUD+6@}mnG(xC^QJjpVyGV-?GoCuX}Ijk!2@JI3Y zCYwsTLv$-MI!31*7J75%qDC&?PHFm^v5?0HP{jJWI4m!hj#&#no68FBu4 zoL=VS;kL`w0gLBzj5bw6dTvy}2zM2E++bQ3OCOGO9r51UN&oX^_AiQ6Q;*~it^K?r zswx`!`|k5cmgTN*y(2z?)EJ2#?j`UZXo1yuLD9VPM#ntqx7-jZ^BSSrvL&Im^gqiE z%>TVK-1F7pmIJHJliZNm7wSKD3sxS~N}YDxp;{EWk!Am7LUO(Exp{bpjUvMhLeU?> z2lpg6TD-rPf7iPxeMT(A{VMyiz{^CHH1;ipqhmBmZ{;|?(Qm=bX;oFue0XwL9R4{Z z(n^e$zmWgdtF#pFUaMhXF_a$u(yL$M%vHy`r0z`pF3sKA`l6HlMpSk4Y2y;*A*r)O z%_MD=mB7@C4&SMJ`L>Q`vrcTNu&cc76!T@t>xAv1q>3PxY{lL8Jdy@a&-Vyz zt{A!KzZPjV?%rr*UR^BU`AgQOKfSoUByVo*K4VtO53#xMrOSJNkg7RbTBX`Og!f7# zOZC>ZI6FVR*RK`BGrJvh>2Pjxc7d7QDd=pvcVH&QCw$2fJ5j{;_kshH4$6;Q zj`q`B@^O|u)9&)kXTN_L=364v_=xky@?hEoeeT0eSNzVIqB7TIoQ7RpdND4y>Yl2J z-t-uK?vUaAmOQT(>FKF6m!{%&+^>r(NgZ!8VQgP^KgS!zYZy&>7G&CYW|_KPBB$A4 zt3WYxuAJ0At@oVtFgkLZ+pMdn{>8w@+^J~Ov!~iU)JM0IM^u$td?9>i4)9Pf*pLF~6`{tQpdDp#Ierqm%RV(}ph{H}|dXXntFJ9@2XE zt?fO*MuCf0MoRs1MPK~oT=P-xZL*SFlxp0sG^2%&j~(R0Vy$Uny?Hn7oIJGTg>lz~W#%-xU$Y8SL zxrlpP?NAMlpRQ^)I<23#1Wt|` zFl*)yHBT%sXndHoI?mZjy6eB@6rN>`w|V<>=W5eV-(Ryzcv+2eOx#_IC56Rut9adq zXJ7rYx)alX(ErA)9^b+72-jZ)19mS*%qzfFPkL04boagJ^ylT9O+k;moA?8Fi+Pvk za~k(#hyNhW=K1uU-H~{3E@CyfbnUM^84Y>UBa@2W^Ay}=9MU{bvku746+Zh>qQr2m z^lQ)fwp&aliW-HQt#ucbN(VIKF3IhRy;3NH1|u@mz4a*zt5UT@yz1S&K{bk z{=7T?V!5WWEcs7^(r)RsolA!wJF0URrQD9%r#xN%msuNI+(dDk^yxM8HD1N5o<78! z-pgM*MQLh3Di~gS@?6NmL+~=YApK-jn1kzS^S7-x{=ShgKQl=mSNh8Z&6c&<*V(F1 zD4xc9?I&NyNijWgFVnnlHlf`iRwp|AIfeCDE)TtvNh{0v=2!Z*3vrR>@X?=(Oxa2~ zzEbH&xxTQkJz^bI@J3KH?wR?k_P*CD$q)QDT<{r79k*k1Yz~#Yiqxj_SBg6ESFvR3 zr!}If9T->-ZhBqHEdHx`V=4Qs zX`teo8lIVGcQidW6fpIr{ej{2Suf_bA8%Kkd4e-l>uoZgZT+73hK--6dEi3-i|qPG z70w%X#8_`W@g4Ye#ac;YLTa+#Za|xZsqr`r*6MBiY3E1le9Q^zgFScr_h^cUS=@0i zLTkis81@^9sg95I|2i;`#QjU70Kc8#f%@3;{07d#XI@p-`oq0F-_i=%BCnMm3NggH z^Bb9^^p!|CgeB%q&YpaCVh8`W(E9y+P~a>t-o57v!;Dh{DSMuAUN1EG(W4r}Myef& z^@Q7%cAOfo+JC;A@UJvyYNT4Bf8(my*teVhwDwlE(i7*C@*mXO(r(m>#znpm6EOV* zHCx4o7+ z7N2cSvOl@}*2envmlQUIXl4duI{XYzVu0FBysJTw|xI6 zsR{dBZ}S;Z7MHt!C{+3Di-kxuvqxH*`MZA*J4UQ}EfUEXZD(_{Uf^l+#2>BU?dMp= z8{L)5uZG^Ne{9wrxg4=c`N-Mvht9(zUVMhx8D)xt%!QIaRV47)rc5@c<~gRK&g2W5 z8GlM2rp@A-K3_A@@=&E}>O8R-q6Pog?5f=Q^j=$CiZ%F_w?HTBj#269WWkwa%@7kw z!xu&M?n>S=x9IFcNaFM7t8TrH`oYed&*LL6B{{&n_FKE>D(9-;>L=05KF>wFBb@z@ znj3{M?&yDUHM*d?(+WdSer=(`Gj0aQbet?eW8kQfJ?V_Fk4Z zU^$|&BluUAgOcfl)OkM5$?*q3E{7Q-%1jx z6+iK-@;lz2H+8gA$F6GCXpnw+bi^4ixzJkPCJ{LBDQ|6YOdxncGY$XE!gIBp?zgtRGzHsGxrv%>$2Aovu>Q?W><+}$kaI@8r6GoQ%#EJk*MJ8si@o> zTlcrWosrHe>iqg?GtYipc&+GUp_6`AO7BXo)VblS;|=aFBhPL%rL|b@6I#XBcdky& z>oMprbdG+1AKrgryNZM90Liu7NbwYFa!?QBtR-D8&lBtP3w2NCD`vm5GJWp&?K-E~ zs~YduhsEf|WBV7%iud15==8TH5;=l{MGx&1xW6B0HDhq}36mmqP~ALmMqzgi-&>{J zBirHr(OFxB&t3IxCC^CE6y<9C>w1!?e83GU%13*+6ujr>qe`lm{yq@8#v2!D)@Kn5 z`k;s{D`;wqigZ@Ixn(0D=t3|XZ0pYIQY!sDujDUIO71u|Sl-;ZaYmtC$oBbEV5tOK zq=rIzO#8O_!>(6VtZ~AW+$pc(>NY~D>f3!|OWpmQQXi`RTQ1)H`NER-t`#+Jsb6n^ z?7^(x-&mzSVWCu}8NJjpAaC=4uUYj_>`#^Ca(}PW|Kdz!b12OkFT8(7sZL^e=fPiD zrFKf~Q=)~sAQOTz(@l2&D1^w~-> zBgZbrtLv-J-NnzBhweW6+I{*DM(*D&_7ydHv`|{Nu0N`=~`}mGNQ34(7qEwZ3B%~t^j^@6- z_RJu>bS<%|I;7WKtEcZ-v5%;qZ9RjHu@}p=4iDIvAC)MZ{`gIKBiY#qDS{g3<5!`e zqoDZgDcF;h(4Jh7)K{WDY;0wGY3jIP`}pggsMT|k;XM-#6J^}ubpc{l8mEQggrD-7 z%Ri@*577r%mq`8H!GFW^v3peF+lneZg$xBA#?reg%C2T@r`zP;R#6^NaCpct?vy*; zo79&5*T-v>Z%e;O$N7EL>N%0cAL{zHFugL?IjJSQk~TK(*R%C6*CM^IeLfnVG+ND! zcd9+1$=RnI(Q(1WPEbg%{>0t(S2`fvqsK~Z)O#BKPkUef4`my^U5}?cDj_Lr>Oq$5 zWH*w^mZnhFK_+TKB7?yYiX@aZS*EgYBiV;WWtosAV;x(TA!ZO`W{mZ{rsw;8|A6=B z_w(@!Gv=P_ysztAj^jA*(-Bvq`#i4SI+0S7el$JA{;w3W;`C{gi?AHyamI(0H}S>@ zeBbAyJ1=JR2WwouhwHa|g60U!2U5So(a6mfwr*T8f5P3HO^%2ypJJ z<~}Y}7DtY^)K;zG>|6eEQ(nA3MdBrN$Zc~>(B}o^w9;|}TDVhyohdhsS_tM&u#vSz z;(e;S=RiRaNbTjW)jGU=y>vpa*t>bO?7}TfC_j?`HQ2F6lPWWWY%rT63NAM>gC1gp zU=e*lEBdcaw~hdKKyPmkl%Lt4g14dl2GH4kIn&WOz|0mymQc*`fO(XS$l518@zdNf zEc1^JT;rMwhz>BIh`BiG>A%HT8u=8=EcOm&^#L%U*n52Gmdo9MvZ>c)nW52c!~JF>JGX%@iAC_)uuHAKzHEkjbH%TK*8h(YzR7!-hXBSH72`9FS1F* zw|FTDNyN)zK0=qvYA1>z=Oy;^0+t==UOR0%`SP?ia_2y;{z}u{Zpb)~IAqL*aaMdP zka>BO+bv%484t=bz9#-z9hvZ^?uOIybNWkzx#{X#(WoL#Bp>ZJ&bQqbz(B4J8kg}s zFIU1$b$+=vb%LQe4(Ls8l()QY^}vq7-Zl+2>U7abln9^oO88}UC6*N|yNs0KZ)H)y znPk{tt$Ac_V-zjjll&e6`15ePnosW+QTH^h}PKfKR58UdH^ggt%br>#LvoX8E2 zw81(S?cxf`a_XKS$C|>Lh#^;vq|KKpcE&nZbYGt{kKI&WR0&bf0tI~T$n+~=v*y%}WDTcZ?95U>WbD~oY@}|73^%p2m@H#w* zE*^7j4jWSc2#x*%{7P7?+j|SuITbfjpp9mvFHLnnITp*&l_sFuU`F0s z!Ea$VtM0yni**hgC3Ej?HSINQRuxn!*bbc^Gam9v|9H92(^N<1o35qs2vYe7Idb^m z7b4#Lo2IS*#R&oLLl?jE%WmZdBDwZhJz>S#%=LDDc&Po|7g}vie3}$$;;q*IY|8_a zy;?4BCgV4oej9fVam2?V##mmf)OGKJg%2Ki{y`kNv9=3nxtQnmdJDH^+ml%5ls4oM z!L1#lSydI5!^d57+#C+*Lt>&2_HGyX&)UJnNARVu1Cr$$D%3n*g<#d86X_L)VCR0q z#aAv5zgFg66@|6>jMPxLR!b#$~}rsAtcHTLup-_^9zJ`M$vOq zfaSpK#r2-gnz4_K%%NnHh~)E*z?|foEHAa-TT?uBkI_z17_y?t>S3>B|n0U=Er z`upMI?}WRuyG1cuI%a2 z@+IyU>Yn{3D*IvfSKSM-q~pqK#xNeKTc^#%POo}<8cjx+IPx5Lj15@Al3o*nM|(|@ zPRF$lkat&9W8fsLQq-gF9JvSdCsSiqK=M zPWYp81j5^j`yZxU#9S!>kBdV!dJXX_!f@*s?@Dn$7Y89}J7fbOt&YJqQ@uy~Sh8OK z3Xj(-5-Bg+li-j~&xZF~nnju*q*OB5nx!a6)10 z(r1XP_Iq^`sa>yEAB~H=JYgKJ8dJg(e;ehEBvx_%a(=sPu45WxDs0E>{|1>R%Cr;B zI$^Wt9ZPtfwLcw#x2U9YpFm{1%GEha{qfZ`>QcDtKYB6JD}_JN$&EbO5yp)U*aala zVsmVX&sZ@{)m`}P^k?pDoql7-9^5%5Tx_KxSvVp}u=Mz0-^>Pz@LMTFNC{S%O*fT7 z24-$-5k zOBpO96>vYq{>Bk)oKS96hU%PC+sPVY&|sIJuyLem)(h!^{N94lD%1AoX8(^AM*g!#7W&F;~G)`Yd24cnj))+%&6mmj=Deu9L zohlMh?MltHFg^v{Nea5hg}#Ck9+5W4U54UQV$qfrrJZ`E1(K-&F%Gg#lw);8Vb5%0 zuoCe`>)NI|pZ#+WEV!!H39^_6g){nXHIhD?ZatATENx@j0$dAVuzKCH*YO=;-a3Rd{Q{eqNBRk! zPY<6g2IGXmO>3cluI5WXoszk`v(~9JAgW}~2+3Y#X*1aI%PH7mUTXmXQJ~-7Wq@HM z_n)7vXbd@peRb9gv9UlbVD3y*H3~ZS)t^}(iT)wz(-=fw#3NhKB=^#t;nrcFGnn`# zW7IFAf8}>ngox$?YDQE4ZpHH*ewxbKWTIP5bRd~Fn3V)}qDqZa4CQ_XWRZWUZtqr> zyHfBq^FH1+zZmObDt_NS0T0tlhSt}m^Kbv%k-nZj5Q7MHz#kvwEoz=aF21hStX3QQ z6dVAaf6gY6>^_C)Kf$&E|6x42=`TP38lX3T4BLOE>B<{x;$4k$wqmQ4VxXwNzW?db7+AQ{^o?rS}Kg2SlZccFgJ50vQmWN0V z*irVyJZ3dPbg({wzK-3?~l+f<(3N*H92_*~gCu)eZtmFOWKQSj9WE}A!_4MUS zYiij}Hqp9g2tl8$9D5#IuGw7A*m(2w@^qU9p~2aVBHTO-T-M9I){n;meoZu@9x*`D z+H0v6yf)iT-QC&hF=}k^v?`q*O90{UdEKcxjcf>7C+sKdx>hFZG?VQaTd(H)L z-fkAfzscUbX{L5Uz}euX!MkYnDhjOY8)f(z)Gs4KWm=H-!6Nh!X_%l6so(6&Ggqn_ z|NW$1Mc#7P`Oco7U5_nrLfOR?Am)**D;@+-eSJDI-Vop%2&1bOR?eX%XDd*48jvpq z%OE;bH8+-9c>L}=+YLonwThgtgGAREM(H?B$X*cs^7B;7&_llHvtZnZcri3M-_~DHA|8bju+W5A zjK3x%i$tvX%aPta7^_>oI4b1D3-^;lgR%f5IzGHSuJ}*iDgdv#2Sbzibr;&IKzRI^ zis1K9y#+#K#|cg8`dM8Z6BMs_ZS?eRT+ zFb)4sv}R{(hqy^9Iqz>zZLH4F_KREhS^iu(V(;>hz^-=08$a($wVe5jNNCVlJfwK1 zMcr}sp5QN9OH%)BL${PG!P8KJ^s3ui?$x+RFPWCAn_4qg4=xx@vgHoWid)Z+AFkn> z);i89so(EC&5zmqs0e5_B#Yv_XOAW|wQ zV(uEuMTSn-OCV@1Q53Dfhbe`UbCFi2Vc*1OeqQer>C1Z88nt8DsiqmA%07WUQErz;PJbJS3aX;%c+fe+nMo5d!bbqrfDY%#w>DKh)o&E<`SO${? zMP(&Q9milZ*FvUM&ZS69Ps#+ZYlnIsX8A6IJA|eqLNXq5D|r+y38}Gz3_q>y^<=4~ z{80X;_jm8nj_PNuZPgt)BO#smDneUmmJ)vB{bs-7ptPjyJV;fZ9(mrjzVR9*w|<1` z;&SQGk!_dWi=wCb)wB<~j$dqxSX%XX2SUKHV*}=U!8Ryv%-?lz>_Jh_-55mpe_fTd|KI zlRrl|z>7=3Ry>;+d~qh1?6jm;nwIMn^!29_BTMddzx!s@GK#I3VY`e++-#=+c}7 z^`*ml$V3h8Q7Ien(@SgU*H16-XfRZS#7K{~B=UH=@VYqY2p61Libt&C;`od^ihZ1s z3;m}R^#rsB@(Nr^;liCXRD;lbsOaMkZBnh{1=f>{MBFStCS=VNh$YC=!~TEABRwh4 z%31d=23tUSpR7SglPaRXk$LEx=ex+;cC}IpL2Rk2H5ZP#?8IdXIYCn;Nl&{ZE3&-4 zZ!N&n-vAef$u<#^SE?PDa;No)T$Yptekt|J)66gDjyTD1yCA;8cU|Ee9GBM2ggd>1 z$Tz?4d)@fh=t}#~M8(2qYIOzuzj#;s-{8LJNU3jEQXvURQjjB5!RgiKo%ZxqC(UjV zMnKZ=ul|!;G3dcOewlI9#C*srbuXXmUAzKq8;e46v$Hm)-zA9q1|i>Ugq{uOok=

!5ST$kuF1Pxfhv+L$kFD6(#FzxEF*vaAHb_f>%tu@c zH$7Jh8&>ba+MdJ5oJ(^Mz2^uC=u?k?lb2dqfFGim$Dt#*(Fl7||vRZ=+=O z2r%b4n-5%AbaY5;NH|MoB9<9Bj%Y&ta(gHJ)H@A#Zq@=Q zk0xrjEkiG#N`f})&6*?G3J^50;KJr`vbd_V*`(!ZgYxE&YY48EYHpt%T3-!onxX7% zGuy>#%&kU7K_&z}91;Wd;GqM#9**taw1saXqb~#IXuNOMmfRSlyx+=c8@%6oMwP{f zyccjJ0QxQ1LQhg-L((mV#{QkB_4o8tlmq&7Rk&1_P|pC6F>+!5g4B?2X7=av1D6qK z)LZwmq6xx1%JUQBLJlY~EUV^@6n@`qtD|zNNC}!}N7Y9Y7e|Nf3Cusc_ zo;`W^|BO_1ZZo@?^LurtWAyIZ#8}XSxNC`lm`fWp2is+MuQ1dGF(QRF8IY zPp1k(mF!WsQj?6mOft>E`4?k-#3Z_$+?!1Mu7j-Z`bB56Ospc$U6Li_(l7h7??d>A zyL#&68#ju|A~&k(JU8e}K+B5EZ@f_Dqdw z(O;dDi#jd$8l{f1)fBY~FrZB=_~w48umB zpo|Hz9(4%#w7n z^XM<(A(cL=8zV!K`wt=VjdP}|O7n)(ze_oaSA((&VpZZiI~)MX@K^YDaS{5oT(vM*;^<@6&`oq~MSr~IknNG- zN|fY|IXY@vfA+Vg+&$|S+98R!GL%&O_4|+GvIMPFGvr3|jy@FDHFuR(^NAb1AX+qb ztk&}B>kF5wxw)s1DX>AqnTs$Et?;IS@$vK*BM3b160`>x}XVvv2`u9^5KoiAp_ zobv(H4A+^np7+XsC|)yF=(V4s*EFw!s2;YQsHB-)GC8xHyg@o0-h-n4XM_M7mQnz0EHN zdKz6yyeqCDq#AE3KdEtG&3O7a@%q{`SzG$99Cs8#WXy80?%xFR%8doR+q~C@ynd3- z>EA9o+NsD*y?3~!z^kF+zhg3LtpUvs;#T^o&@bDNo0o^X3r<#+d%4_FQtr(TqSVDK z-2Ruy@J*Ek-`!e>+k4B5GvOLf@PZqrp;q=zX&17pOfy6Iz$FG{Lk#{t#a&v5JsuYU?&<(N0` z4rk4kKql%dYV($_CeP3bj=%;dXRxtX7jCpv#L7A~$38gBaqmsXq3v6$dFsI_2#jQz8JiOF>#|v6 z)HarVfTkB><4x>wg*AJ`cpXbSK`t+UAqB+}h`v1h`}*Y?BMtV8t1T}_{t zN0W`BkKV*E0)HV z=s@iV72tN2<7O`OvI!o8ZDr8fv_Tzm4lCPL&4A%UP>e5};y}!+JPsO6Z=pwH380VA zU<+Kq{iKppCaAN@H_4aTBG&ObAGklLy%v?mKU`yU0kxvjPM%Pw^aLxk$f{N!#@2>1 z+C;pqC*8s2=q9$`fW6}R z;B2r(tz1+Jju+F_*czSO)z~v2FjROFuWGs&fITOR; zBJZ*JlOCp;?G;>ps^OUJZklsweR=j@gR4z@TTL|Ndc8)P_R7eWIg2{%jP}qSpa;X% zNR=_cbmJa9{R4!bU=S^FTU0KE??fDh7r{{7vPg`&h)Z6xTzibuagASp<+lcUsYlaz z3e;75%q@cs^?ifTwL1=1X`aok?C08=ZQ?_ytez>~TP29k(|bix$57KT{L(rFKVqrC z!MCUwLjdSrdTwf%P1n%2WsTEiTEpqBlTu=i(W1jSk?f|0T;%(;bm1U#v0`iLE*RI_ zEkZ0(Pn}KCTQ4dFat7{9^$VIvRsJx!T$SCtD^KyZ?j!?LV(F%gagpbpw^Yn!$LuaV7S!pdoPh~?7LS2=8DcQbjYJi@D+8tRHc z5cvK1yOY3!x*M2~8lq&B$O}ZUpLLxh`;x1;*-aUvalZS=;ZJZHnuR8b`Cr2IbzZXD zJ1i?}{UWrTjt$6{#j^F%I(2pJvnc>i3~l~LhV%NVmy#CEY~lES z>QrQlM~tr?qSKJw<3UnC^1F+Z{xA{~86=6Y2X^FWkJ9+_#;U^g;?H9a-k)xw*(zcU zW5j3N|GZ@ZXfi83hOk}DpaHhH$&CHFWcAj@OhO6t&%V2U@phXHjCEUjl20&b0RRuK1mY6Hl8=5na2RQiw#@^VppG%*=yO>;NZGgSsvqUywe%V94k*Q<3O**R>ddb5TZWXD-*-lr}`8o7y1 z1>bQ_@bPcZRS&+790R1)MR?QEXd3!gxaO@-k_g>i)|bZ?qWVGoHUd}2z_RB9a(JE^ z!bVj;#L{$-%d3lpuhj-#`5S2B7~n&HsZXGDoO;`%heGA+Mz9 SeP_3ZnO?qj3475c;{N~?2;#5+ literal 0 HcmV?d00001 diff --git a/docs/src/model_docs/lateral/kinwave.md b/docs/src/model_docs/lateral/kinwave.md index 3e7b92727..1fdde084a 100644 --- a/docs/src/model_docs/lateral/kinwave.md +++ b/docs/src/model_docs/lateral/kinwave.md @@ -55,30 +55,64 @@ kinematic wave for surface water routing, as a cyclic parameter or as part of fo also [Input section](@ref)). ## Subsurface flow routing -In the SBM model the kinematic wave approach is used to route subsurface flow laterally. The -saturated store ``S`` can be drained laterally by saturated downslope subsurface flow per -unit width of slope ``w`` [m] according to: +In the SBM model the kinematic wave approach is used to route subsurface flow laterally. +Different vertical hydraulic conductivity depth profiles are possible as part of the +vertical [SBM](@ref soil) concept, and these profiles (after unit conversion) are also used +to compute lateral subsurface flow. The following profiles (see [SBM](@ref soil) for a +detailed description) are available: +- `exponential` (default) +- `exponential_constant` +- `layered` +- `layered_exponential` +For the profiles `exponential` and `exponential_constant`, the saturated store ``S`` is +drained laterally by saturated downslope subsurface flow for a slope with width ``w`` [m] +according to: ```math - q=\frac{K_{0}\mathit{tan(\beta)}}{f}(e^{(-fz_{i})}-e^{(-fz_{t})}) + Q = \begin{cases} + \frac{K_0\tan(\beta)}{f}\left(e^{(-fz_{i})}-e^{(-fz_\mathrm{exp})}\right) w + + K_0e^{(-fz_\mathrm{exp})}(z_t-z_\mathrm{exp})\tan(\beta) w & \text{if $z_i < z_\mathrm{exp}$}\\ + \\ + K_0e^{(-fz_\mathrm{exp})}(z_t - z_i)\tan(\beta) w & \text{if $z_i \ge z_\mathrm{exp}$}, + \end{cases} ``` -where ``\beta`` is element slope angle [deg.], ``q`` is subsurface flow [m``^{2}``/t], -``K_{0}`` is the saturated hydraulic conductivity at the soil surface [m/t], ``z_{i}`` is -the water table depth [m], ``z_{t}`` is total soil depth [m], and ``f`` is a scaling -parameter [m``^{-1}``], that controls the decrease of vertical saturated conductivity with -depth. +where ``\beta`` is element slope angle, ``Q`` is subsurface flow [m``^{3}`` d``^{-1}``], +``K_0`` is the saturated hydraulic conductivity at the soil surface [m d``^{-1}``], ``z_i`` +is the water table depth [m], ``z_{t}`` is the total soil depth [m], ``f`` is a scaling +parameter [m``^{-1}``] that controls the decrease of ``K_0`` with depth and +``z_\mathrm{exp}`` [m] is the depth from soil surface for which the exponential decline of +``K_0`` is valid. For the `exponential` profile, ``z_\mathrm{exp}`` is equal to ``z_t``. Combining with the following continuity equation: ```math - (\theta_s-\theta_r)\frac{\partial h}{\partial t} = -w\frac{\partial q}{\partial x} + wr + (\theta_s-\theta_r)w\frac{\partial h}{\partial t} = -\frac{\partial Q}{\partial x} + wr ``` -where ``h`` is the water table height [m], ``x`` is the distance downslope [m], and ``r`` -is the net input rate [m/t] to the saturated store. Substituting for ``h (\frac{\partial -q}{\partial h})``, gives: +where ``h`` is the water table height [m], ``x`` is the distance downslope [m], and ``r`` is +the net input rate [m d``^{-1}``] to the saturated store. Substituting for ``h +(\frac{\partial Q}{\partial h})``, gives: ```math - w \frac{\partial q}{\partial t} = -cw\frac{\partial q}{\partial x} + cwr + \frac{\partial Q}{\partial t} = -c\frac{\partial Q}{\partial x} + cwr ``` -where celerity ``c = \frac{K_{0}\mathit{tan(\beta)}}{(\theta_s-\theta_r)} e^{(-fz_{i})}`` +where celerity ``c`` is calculated as follows: +```math + c = \begin{cases} + \frac{K_0e^{(-fz_{i})}\tan(\beta)}{(\theta_s-\theta_r)} + + \frac{K_0e^{(-fz_\mathrm{exp})}\tan(\beta)}{(\theta_s-\theta_r)} & \text{if $z_i < z_\mathrm{exp}$}\\ + \\ + \frac{K_0e^{(-fz_\mathrm{exp})}\tan(\beta)}{(\theta_s-\theta_r)} & \text{if $z_i \ge z_\mathrm{exp}$}. + \end{cases} +``` + +For the `layered` and `layered_exponential` profiles the equivalent horizontal hydraulic +conductivity ``K_h`` [m d``^{-1}``] is calculated for water table height ``h = z_t-z_i`` +[m], and lateral subsurface flow is calculated as follows: +```math + Q = K_h h \tan(\beta) w, +``` +and celerity ``c`` is given by: +```math + c = \frac{K_h \tan(\beta)}{(\theta_s-\theta_r)}. +``` The kinematic wave equation for lateral subsurface flow is solved iteratively using Newton's method. diff --git a/docs/src/model_docs/params_lateral.md b/docs/src/model_docs/params_lateral.md index 97234d044..7ad19afaa 100644 --- a/docs/src/model_docs/params_lateral.md +++ b/docs/src/model_docs/params_lateral.md @@ -155,11 +155,11 @@ between parentheses. The Table below shows the parameters (fields) of struct `LateralSSF`, including a description of these parameters, the unit, and default value if applicable. The parameters in bold represent model parameters that can be set through static input data (netCDF). The -soil related parameters `f`, `soilthickness`, `θₛ` and `θᵣ` are derived from the vertical -`SBM` concept (including unit conversion for `f` and `soilthickness`), and can be listed in -the TOML configuration file under `[input.vertical]`, to map the internal model parameter to -the external netCDF variable. The internal slope model parameter `βₗ` is set through the -TOML file as follows: +soil related parameters `f`, `soilthickness`, `z_exp`, `θₛ` and `θᵣ` are derived from the +vertical `SBM` concept (including unit conversion for `f`, `z_exp` and `soilthickness`), and +can be listed in the TOML configuration file under `[input.vertical]`, to map the internal +model parameter to the external netCDF variable. The internal slope model parameter `βₗ` is +set through the TOML file as follows: ```toml [input.lateral.land] @@ -168,22 +168,33 @@ slope = "Slope" The parameter `kh₀` is computed by multiplying the vertical hydraulic conductivity at the soil surface `kv₀` (including unit conversion) of the vertical `SBM` concept with the -external parameter `ksathorfrac` \[-\] (default value of 1.0). The external parameter -`ksathorfrac` can be set as follows through the TOML file: +internal parameter `khfrac` \[-\] (default value of 1.0). The internal model parameter +`khfrac` is set through the TOML file as follows: ```toml [input.lateral.subsurface] ksathorfrac = "KsatHorFrac" ``` -The `ksathorfrac` parameter compensates for anisotropy, small scale `kv₀` measurements (soil +The `khfrac` parameter compensates for anisotropy, small scale `kv₀` measurements (soil core) that do not represent larger scale hydraulic conductivity, and smaller flow length scales (hillslope) in reality, not represented by the model resolution. +For the vertical [SBM](@ref params_sbm) concept different vertical hydraulic conductivity +depth profiles are possible, and these also determine which `LateralSSF` parameters are used +including the input requirements for the computation of lateral subsurface flow. For the +`exponential` profile the model parameters `kh₀` and `f` are used. For the +`exponential_constant` profile `kh₀` and `f` are used, and `z_exp` is required as part of +`[input.vertical]`. For the `layered` profile, `SBM` model parameter `kv` is used, and for +the `layered_exponential` profile `kv` is used and `z_exp` is required as part of +`[input.vertical]`. + | parameter | description | unit | default | |:---------------| --------------- | ---------------------- | ----- | | `kh₀` | horizontal hydraulic conductivity at soil surface | m d``^{-1}`` | 3.0 | -| **`f`** | a scaling parameter (controls exponential decline of kh₀) | m``^{-1}`` | 1.0 | +| **`f`** | a scaling parameter (controls exponential decline of `kh₀`) | m``^{-1}`` | 1.0 | +| `kh` | horizontal hydraulic conductivity | m d``^{-1}`` | - | +| **`khfrac`** (`ksathorfrac`) | a muliplication factor applied to vertical hydraulic conductivity `kv` | - | 100.0 | | **`soilthickness`** | soil thickness | m | 2.0 | | **`θₛ`** | saturated water content (porosity) | - | 0.6 | | **`θᵣ`** | residual water content | - | 0.01 | @@ -192,13 +203,13 @@ scales (hillslope) in reality, not represented by the model resolution. | `dl` | drain length | m | - | | `dw` | drain width | m | - | | `zi` | pseudo-water table depth (top of the saturated zone) | m | - | +| **`z_exp`** | depth from soil surface for which exponential decline of `kh₀` is valid | m | - | | `exfiltwater` | exfiltration (groundwater above surface level, saturated excess conditions) | m Δt⁻¹ | - | | `recharge` | net recharge to saturated store | m``^2`` Δt⁻¹ | - | | `ssf` | subsurface flow | m``^3`` d``{-1}`` | - | | `ssfin` | inflow from upstream cells | m``^3`` d``{-1}`` | - | | `ssfmax` | maximum subsurface flow | m``^2`` d``{-1}`` | - | | `to_river` | part of subsurface flow that flows to the river | m``^3`` d``{-1}`` | - | -| `wb_pit` | boolean location (0 or 1) of a waterbody (wb, reservoir or lake) | - | - | ## Local inertial diff --git a/docs/src/model_docs/params_vertical.md b/docs/src/model_docs/params_vertical.md index 03a7817c0..fde92c15b 100644 --- a/docs/src/model_docs/params_vertical.md +++ b/docs/src/model_docs/params_vertical.md @@ -16,6 +16,20 @@ external netCDF variable `Sl`: specific_leaf = "Sl" ``` +Different [vertical hydraulic conductivity depth profiles](@ref soil): `exponential` +(default), `exponential_constant`, `layered` and `layered_exponential` can be provided +through the TOML file. Below an example for the `exponential_constant` profile: + +```toml +[input.vertical] +ksat_profile = "exponential_constant" +``` + +For the `exponential` profile the input parameters `kv₀` and `f` are used. For the +`exponential_constant` profile `kv₀` and `f` are used, and `z_exp` is required as input. For +the `layered` profile, input parameter `kv` is used, and for the `layered_exponential` +profile `kv` is used and `z_layered` is required as input. + | parameter | description | unit | default | |:---------------| --------------- | ---------------------- | ----- | | **`cfmax`** | degree-day factor | mm ᵒC``^{-1}`` Δt``^{-1}`` | 3.75653 mm ᵒC``^{-1}`` day``^{-1}`` | @@ -33,7 +47,10 @@ specific_leaf = "Sl" | **`θₛ`** (`theta_s`) | saturated water content (porosity) | - | 0.6 | | **`θᵣ`** (`theta_r`) | residual water content | - | 0.01 | | **`kv₀`** (`kv_0`) | Vertical hydraulic conductivity at soil surface | mm Δt``^{-1}`` | 3000.0 mm day``^{-1}``| -| **`f`** | scaling parameter (controls exponential decline of kv₀) | mm``^{-1}`` | 0.001 | +| **`kv`** | Vertical hydraulic conductivity per soil layer | mm Δt``^{-1}`` | 1000.0 mm day``^{-1}``| +| **`f`** | scaling parameter (controls exponential decline of `kv₀`) | mm``^{-1}`` | 0.001 | +| **`z_exp`** | Depth from soil surface for which exponential decline of `kv₀` is valid | mm | - | +| **`z_layered`** | Depth from soil surface for which layered profile (of `layered_exponential`) is valid | mm | - | | **`hb`** | air entry pressure of soil (Brooks-Corey) | cm | 10.0 | | **`soilthickness`** | soil thickness | mm | 2000.0 | | **`infiltcappath`** | infiltration capacity of the compacted areas | mm Δt``^{-1}`` | 10.0 mm day``^{-1}`` | @@ -59,6 +76,7 @@ specific_leaf = "Sl" | `n` | number of grid cells | - | - | | `nlayers` | number of soil layers | - | - | | `n_unsatlayers` | number of unsaturated soil layers | - | - | +| `nlayers_kv` | number of soil layers with vertical hydraulic conductivity value `kv` | - | - | | `riverfrac` | fraction of river | - | - | | `act_thickl` | thickness of soil layers | mm | - | | `sumlayers` | cumulative sum of soil layers thickness, starting at soil surface | mm | - | @@ -69,7 +87,6 @@ specific_leaf = "Sl" | `zi`| pseudo-water table depth (top of the saturated zone) | mm | - | | `soilwatercapacity`| soilwater capacity | mm | - | | `canopystorage`| canopy storage | mm | - | -| `canopygapfraction` | canopygapfraction | - | - | |**`precipitation`** | precipitation | mm Δt``^{-1}``| - | | **`temperature`** | temperature | ᵒC | - | | **`potential_evaporation`** | potential evaporation | mm Δt``^{-1}`` | - | @@ -111,7 +128,6 @@ specific_leaf = "Sl" | `snow` | snow storage | mm | - | | `snowwater` | liquid water content in the snow pack | mm | - | | `rainfallplusmelt` | snowmelt + precipitation as rainfall | mm Δt``^{-1}`` | - | -| `glacierstore` | water within the glacier | mm | - | | `tsoil` | top soil temperature | ᵒC | - | | **`leaf_area_index`** | leaf area index | m``^2`` m``{-2}`` | - | | `waterlevel_land` | water level land | mm | - | diff --git a/docs/src/model_docs/vertical/sbm.md b/docs/src/model_docs/vertical/sbm.md index 14b11338d..7f174a441 100644 --- a/docs/src/model_docs/vertical/sbm.md +++ b/docs/src/model_docs/vertical/sbm.md @@ -350,14 +350,28 @@ t``^{-1}``]) is in that case controlled by the saturated hydraulic conductivity ```math st=K_{\mathit{sat}}\frac{U_{s}}{S_{d}} ``` -Saturated conductivity (``K_{sat}`` [mm t``^{-1}``]) declines with soil depth (``z`` [mm]) -in the model according to: + +Four different saturated hydraulic conductivity depth profiles (`ksat_profile`) are +available and a `ksat_profile` can be specified in the TOML file as follows: + +```toml +[input.vertical] +ksat_profile = "exponential_constant" # optional, one of ("exponential", "exponential_constant", "layered", "layered_exponential"), default is "exponential" +``` + +Soil measurements are often available for about the upper 1.5-2 m of the soil column to +estimate the saturated hydraulic conductivity, while these measurements are often lacking +for soil depths beyond 1.5-2 m. These different profiles allow to extent the saturated +hydraulic conductivity profile based on measurements (either an exponential fit or hydraulic +conductivity value per soil layer) with an exponential or constant profile. By default, with +`ksat_profile` "exponential", the saturated hydraulic conductivity (``K_{sat}`` [mm +t``^{-1}``]) declines with soil depth (``z`` [mm]) in the model according to: ```math - K_{sat}=K_{0}e^{(-fz)} + K_{sat}=K_{0}e^{(-fz)}, ``` -where ``K_{0}`` [mm t``^{-1}``] is the saturated conductivity at the soil surface and ``f`` -is a scaling parameter [mm``^{-1}``]. +where ``K_{0}`` [mm t``^{-1}``] is the saturated hydraulic conductivity at the soil surface +and ``f`` is a scaling parameter [mm``^{-1}``]. The plot below shows the relation between soil depth ``z`` and saturated hydraulic conductivity ``K_{sat}`` for different values of ``f``. @@ -385,6 +399,30 @@ conductivity ``K_{sat}`` for different values of ``f``. end # hide ``` +With `ksat_profile` "exponential\_constant", ``K_{sat}`` declines exponentially with soil +depth ``z`` until ``z_\mathrm{exp}`` [mm] below the soil surface, and stays constant at and +beyond soil depth ``z_\mathrm{exp}``: + +```math + K_{sat} = \begin{cases} + K_{0}e^{(-fz)} & \text{if $z < z_\mathrm{exp}$}\\ + K_{0}e^{(-fz_\mathrm{exp})} & \text{if $z \ge z_\mathrm{exp}$}. + \end{cases} +``` + +It is also possible to provide a ``K_{sat}`` value per soil layer by specifying +`ksat_profile` "layered", these ``K_{sat}`` values are used directly to compute the vertical +transfer of water between soil layers and to the saturated store ``S``. Finally, with the +`ksat_profile` "layered\_exponential" a ``K_{sat}`` value per soil layer is used until depth +``z_\mathrm{layered}`` below the soil surface, and beyond ``z_\mathrm{layered}`` an +exponential decline of ``K_{sat}`` (of the soil layer with bottom ``z_\mathrm{layered}``) +controlled by ``f`` occurs. The different available `ksat_profle` options are schematized in +the figure below where the blue line represents the ``K_{sat}`` value. + +![ksat_profiles](../../images/sbm_ksat_profiles.png) + +*Overview of available `ksat_profile` options, for a soil column with five layers* + ### Infiltration The water available for infiltration is taken as the rainfall including meltwater. diff --git a/src/flow.jl b/src/flow.jl index c9c8d25e8..d7dd800ac 100644 --- a/src/flow.jl +++ b/src/flow.jl @@ -389,6 +389,8 @@ end @get_units @exchange @grid_type @grid_location @with_kw struct LateralSSF{T} kh₀::Vector{T} | "m d-1" # Horizontal hydraulic conductivity at soil surface [m d⁻¹] f::Vector{T} | "m-1" # A scaling parameter [m⁻¹] (controls exponential decline of kh₀) + kh::Vector{T} | "m d-1" # Horizontal hydraulic conductivity [m d⁻¹] + khfrac::Vector{T} | "-" # A muliplication factor applied to vertical hydraulic conductivity `kv` [-] soilthickness::Vector{T} | "m" # Soil thickness [m] θₛ::Vector{T} | "-" # Saturated water content (porosity) [-] θᵣ::Vector{T} | "-" # Residual water content [-] @@ -397,6 +399,7 @@ end dl::Vector{T} | "m" # Drain length [m] dw::Vector{T} | "m" # Flow width [m] zi::Vector{T} | "m" # Pseudo-water table depth [m] (top of the saturated zone) + z_exp::Vector{T} | "m" # Depth [m] from soil surface for which exponential decline of kv₀ is valid exfiltwater::Vector{T} | "m Δt-1" # Exfiltration [m Δt⁻¹] (groundwater above surface level, saturated excess conditions) recharge::Vector{T} | "m2 Δt-1" # Net recharge to saturated store [m² Δt⁻¹] ssf::Vector{T} | "m3 d-1" # Subsurface flow [m³ d⁻¹] @@ -411,7 +414,7 @@ end end -function update(ssf::LateralSSF, network, frac_toriver) +function update(ssf::LateralSSF, network, frac_toriver, ksat_profile) @unpack subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes = network ns = length(subdomain_order) @@ -433,22 +436,42 @@ function update(ssf::LateralSSF, network, frac_toriver) upstream_nodes[n], eltype(ssf.to_river), ) - - ssf.ssf[v], ssf.zi[v], ssf.exfiltwater[v] = kinematic_wave_ssf( - ssf.ssfin[v], - ssf.ssf[v], - ssf.zi[v], - ssf.recharge[v], - ssf.kh₀[v], - ssf.βₗ[v], - ssf.θₛ[v] - ssf.θᵣ[v], - ssf.f[v], - ssf.soilthickness[v], - ssf.Δt, - ssf.dl[v], - ssf.dw[v], - ssf.ssfmax[v], - ) + if (ksat_profile == "exponential") || + (ksat_profile == "exponential_constant") + ssf.ssf[v], ssf.zi[v], ssf.exfiltwater[v] = kinematic_wave_ssf( + ssf.ssfin[v], + ssf.ssf[v], + ssf.zi[v], + ssf.recharge[v], + ssf.kh₀[v], + ssf.βₗ[v], + ssf.θₛ[v] - ssf.θᵣ[v], + ssf.f[v], + ssf.soilthickness[v], + ssf.Δt, + ssf.dl[v], + ssf.dw[v], + ssf.ssfmax[v], + ssf.z_exp[v], + ksat_profile, + ) + elseif (ksat_profile == "layered") || + (ksat_profile == "layered_exponential") + ssf.ssf[v], ssf.zi[v], ssf.exfiltwater[v] = kinematic_wave_ssf( + ssf.ssfin[v], + ssf.ssf[v], + ssf.zi[v], + ssf.recharge[v], + ssf.kh[v], + ssf.βₗ[v], + ssf.θₛ[v] - ssf.θᵣ[v], + ssf.soilthickness[v], + ssf.Δt, + ssf.dl[v], + ssf.dw[v], + ssf.ssfmax[v], + ) + end end end end diff --git a/src/horizontal_process.jl b/src/horizontal_process.jl index b61585608..380edd39e 100644 --- a/src/horizontal_process.jl +++ b/src/horizontal_process.jl @@ -68,8 +68,65 @@ function kin_wave!(Q, graph, toposort, Qold, q, α, β, DCL, Δt) return Q end -"Kinematic wave for lateral subsurface flow for a single cell and timestep" -function kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh₀, β, θₑ, f, d, Δt, Δx, dw, ssfmax) +"Returns water table depth `zi` based on lateral subsurface flow `ssf` and hydraulic conductivity profile `ksat_profile`" +function ssf_water_table_depth(ssf, kh₀, β, f, d, dw, z_exp, ksat_profile) + if ksat_profile == "exponential" + zi = log((f * ssf) / (dw * kh₀ * β) + exp(-f * d)) / -f + elseif ksat_profile == "exponential_constant" + ssf_constant = kh₀ * β * exp(-f * z_exp) * (d - z_exp) * dw + if ssf > ssf_constant + zi = log((f * (ssf - ssf_constant)) / (dw * kh₀ * β) + exp(-f * z_exp)) / -f + else + zi = d - ssf / (dw * kh₀ * β * exp(-f * z_exp)) + end + end + return zi +end + +"Returns kinematic wave celecity `Cn` of lateral subsurface flow based on hydraulic conductivity profile `ksat_profile`" +function ssf_celerity(zi, kh₀, β, θₑ, f, z_exp, ksat_profile) + if ksat_profile == "exponential" + Cn = (kh₀ * exp(-f * zi) * β) / θₑ + elseif ksat_profile == "exponential_constant" + Cn_const = (kh₀ * exp(-f * z_exp) * β) / θₑ + if zi < z_exp + Cn = (kh₀ * exp(-f * zi) * β) / θₑ + Cn_const + else + Cn = Cn_const + end + end + return Cn +end + +""" + kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh₀, β, θₑ, f, d, Δt, Δx, dw, ssfmax, z_exp, ksat_profile) + +Kinematic wave for lateral subsurface flow for a single cell and timestep. An exponential +decline of hydraulic conductivity at the soil surface `kh₀`, controllled by parameter `f`, +is assumed. The hydraulic conductivity profile `ksat_profile` is either `exponential` or +`exponential_constant`, with `z_exp` the depth from the soil surface for which the +exponential decline of `kh₀` is valid. + +Returns lateral subsurface flow `ssf`, water table depth `zi` and exfiltration rate +`exfilt`. +""" +function kinematic_wave_ssf( + ssfin, + ssfₜ₋₁, + ziₜ₋₁, + r, + kh₀, + β, + θₑ, + f, + d, + Δt, + Δx, + dw, + ssfmax, + z_exp, + ksat_profile, +) ϵ = 1.0e-12 max_iters = 3000 @@ -81,10 +138,10 @@ function kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh₀, β, θ ssf = (ssfₜ₋₁ + ssfin) / 2.0 count = 1 - # Estimate zi on the basis of the relation between subsurfacel flow and zi - zi = log((f * ssf) / (dw * kh₀ * β) + exp(-f * d)) / -f + # Estimate zi on the basis of the relation between subsurface flow and zi + zi = ssf_water_table_depth(ssf, kh₀, β, f, d, dw, z_exp, ksat_profile) # Reciprocal of derivative delta Q/ delta z_i, constrained w.r.t. neff on the basis of the continuity equation) - Cn = (kh₀ * exp(-f * zi) * β) / θₑ + Cn = ssf_celerity(zi, kh₀, β, θₑ, f, z_exp, ksat_profile) # Term of the continuity equation for Newton-Raphson iteration for iteration 1 # because celerity Cn is depending on zi, the increase or decrease of zi is moved to the recharge term of the continuity equation # then (1./Cn)*ssfₜ₋₁ can be replaced with (1./Cn)*ssf, and thus celerity and lateral flow rate ssf are then in line @@ -104,10 +161,9 @@ function kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh₀, β, θ # Start while loop of Newton-Raphson iteration m until continuity equation approaches zero while true # Estimate zi on the basis of the relation between lateral flow rate and groundwater level - zi = log((f * ssf) / (dw * kh₀ * β) + exp(-f * d)) / -f + zi = ssf_water_table_depth(ssf, kh₀, β, f, d, dw, z_exp, ksat_profile) # Reciprocal of derivative delta Q/ delta z_i, constrained w.r.t. neff on the basis of the continuity equation - Cn = (kh₀ * exp(-f * zi) * β) / θₑ - + Cn = ssf_celerity(zi, kh₀, β, θₑ, f, z_exp, ksat_profile) # Term of the continuity equation for given Newton-Raphson iteration m # because celerity Cn is depending on zi, the increase or decrease of zi is moved to the recharge term of the continuity equation # then (1./Cn)*ssfₜ₋₁ can be replaced with (1./Cn)*ssf, and thus celerity and lateral flow rate ssf are then in line @@ -142,6 +198,71 @@ function kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh₀, β, θ end end +""" + kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh, β, θₑ, d, Δt, Δx, dw, ssfmax) + +Kinematic wave for lateral subsurface flow for a single cell and timestep, based on +(average) hydraulic conductivity `kh`. + +Returns lateral subsurface flow `ssf`, water table depth `zi` and exfiltration rate +`exfilt`. +""" +function kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh, β, θₑ, d, Δt, Δx, dw, ssfmax) + + ϵ = 1.0e-12 + max_iters = 3000 + + if ssfin + ssfₜ₋₁ ≈ 0.0 && r <= 0.0 + return 0.0, d, 0.0 + else + # initial estimate + ssf = (ssfₜ₋₁ + ssfin) / 2.0 + count = 1 + # celerity (Cn) + Cn = (β * kh) / θₑ + # constant term of the continuity equation for Newton-Raphson + c = (Δt / Δx) * ssfin + (1.0 / Cn) * ssfₜ₋₁ + r + # continuity equation of which solution should be zero + fQ = (Δt / Δx) * ssf + (1.0 / Cn) * ssf - c + # Derivative of the continuity equation + dfQ = (Δt / Δx) + 1.0 / Cn + # Update lateral subsurface flow estimate ssf + ssf = ssf - (fQ / dfQ) + if isnan(ssf) + ssf = 0.0 + end + ssf = max(ssf, 1.0e-30) + + # Start while loop of Newton-Raphson iteration + while true + fQ = (Δt / Δx) * ssf + (1.0 / Cn) * ssf - c + dfQ = (Δt / Δx) + 1.0 / Cn + ssf = ssf - (fQ / dfQ) + if isnan(ssf) + ssf = 0.0 + end + ssf = max(ssf, 1.0e-30) + if (abs(fQ) <= ϵ) || (count >= max_iters) + break + end + count += 1 + end + + # Constrain the lateral subsurface flow rate ssf + ssf = min(ssf, (ssfmax * dw)) + # On the basis of the lateral flow rate, estimate the amount of groundwater level above surface (saturation excess conditions), then rest = negative + zi = ziₜ₋₁ - (ssfin * Δt + r * Δx - ssf * Δt) / (dw * Δx) / θₑ + if zi > d + ssf = max(ssf - (dw * Δx) * θₑ * (zi - d), 1.0e-30) + end + # Exfiltration rate and set zi to zero. + exfilt = min(zi, 0.0) * -θₑ + zi = clamp(zi, 0.0, d) + + return ssf, zi, exfilt + end +end + """ accucapacitystate!(material, network, capacity) -> material diff --git a/src/sbm.jl b/src/sbm.jl index 77c44347c..adc103df4 100644 --- a/src/sbm.jl +++ b/src/sbm.jl @@ -9,6 +9,8 @@ nlayers::Vector{Int} | "-" # Number of unsaturated soil layers n_unsatlayers::Vector{Int} | "-" + # Number of soil layers with vertical hydraulic conductivity value `kv` + nlayers_kv::Vector{Int} | "-" # Fraction of river [-] riverfrac::Vector{T} | "-" # Saturated water content (porosity) [-] @@ -17,6 +19,8 @@ θᵣ::Vector{T} | "-" # Vertical hydraulic conductivity [mm Δt⁻¹] at soil surface kv₀::Vector{T} + # Vertical hydraulic conductivity [mm Δt⁻¹] per soil layer + kv::Vector{SVector{N,T}} | "-" # Muliplication factor [-] applied to kv_z (vertical flow) kvfrac::Vector{SVector{N,T}} | "-" # Air entry pressure [cm] of soil (Brooks-Corey) @@ -55,6 +59,10 @@ throughfall::Vector{T} # A scaling parameter [mm⁻¹] (controls exponential decline of kv₀) f::Vector{T} | "mm-1" + # Depth [mm] from soil surface for which exponential decline of kv₀ is valid + z_exp::Vector{T} | "mm" + # Depth [mm] from soil surface for which layered profile is valid + z_layered::Vector{T} | "mm" # Amount of water in the unsaturated store, per layer [mm] ustorelayerdepth::Vector{SVector{N,T}} | "mm" # Saturated store [mm] @@ -254,6 +262,7 @@ function initialize_sbm(nc, config, riverfrac, inds) Δt = Second(config.timestepsecs) config_thicknesslayers = get(config.model, "thicknesslayers", Float[]) + ksat_profile = get(config.input.vertical, "ksat_profile", "exponential")::String if length(config_thicknesslayers) > 0 thicknesslayers = SVector(Tuple(push!(Float.(config_thicknesslayers), mv))) sumlayers = pushfirst(cumsum(thicknesslayers), 0.0) @@ -463,9 +472,15 @@ function initialize_sbm(nc, config, riverfrac, inds) cmax, e_r, canopygapfraction, sl, swood, kext = initialize_canopy(nc, config, inds) + θₑ = θₛ .- θᵣ + soilwatercapacity = soilthickness .* θₑ + satwaterdepth = 0.85 .* soilwatercapacity # cold state value for satwaterdepth + zi = max.(0.0, soilthickness .- satwaterdepth ./ θₑ) # cold state value for zi + # these are filled in the loop below # TODO see if we can replace this approach nlayers = zeros(Int, n) + n_unsatlayers = zeros(Int, n) act_thickl = zeros(Float, maxlayers, n) s_layers = zeros(Float, maxlayers + 1, n) @@ -478,39 +493,93 @@ function initialize_sbm(nc, config, riverfrac, inds) nlayers[i] = nlayers_ act_thickl[:, i] = act_thickl_ s_layers[:, i] = s_layers_ + _, n_unsatlayers[i] = set_layerthickness(zi[i], sumlayers, thicknesslayers) else nlayers[i] = 1 act_thickl[:, i] = SVector(soilthickness[i]) s_layers[:, i] = pushfirst(cumsum(SVector(soilthickness[i])), 0.0) end end - - # needed for derived parameters below act_thickl = svectorscopy(act_thickl, Val{maxlayers}()) - θₑ = θₛ .- θᵣ - soilwatercapacity = soilthickness .* θₑ - satwaterdepth = 0.85 .* soilwatercapacity # cold state value for satwaterdepth - zi = max.(0.0, soilthickness .- satwaterdepth ./ θₑ) # cold state value for zi # copied to array of sarray below vwc = fill(mv, maxlayers, n) vwc_perc = fill(mv, maxlayers, n) + sumlayers = svectorscopy(s_layers, Val{maxlayers + 1}()) + + # ksat profiles + if ksat_profile == "exponential" + z_exp = soilthickness + z_layered = fill(mv, n) + kv = fill(mv, (maxlayers, n)) + nlayers_kv = fill(0, n) + elseif ksat_profile == "exponential_constant" + z_exp = + ncread(nc, config, "vertical.z_exp"; optional = false, sel = inds, type = Float) + z_layered = fill(mv, n) + kv = fill(mv, (maxlayers, n)) + nlayers_kv = fill(0, n) + elseif ksat_profile == "layered" || ksat_profile == "layered_exponential" + z_exp = fill(mv, n) + kv = + ncread( + nc, + config, + "vertical.kv"; + sel = inds, + defaults = 1000.0, + type = Float, + dimname = :layer, + ) .* (Δt / basetimestep) + if size(kv, 1) != maxlayers + parname = param(config.input.vertical, "kv") + size1 = size(kv, 1) + error("$parname needs a layer dimension of size $maxlayers, but is $size1") + end + if ksat_profile == "layered" + z_layered = soilthickness + nlayers_kv = nlayers + else + z_layered = ncread( + nc, + config, + "vertical.z_layered"; + optional = false, + sel = inds, + type = Float, + ) + nlayers_kv = fill(0, n) + for i in eachindex(nlayers_kv) + layers = @view sumlayers[i][2:nlayers[i]] + _, k = findmin(abs.(z_layered[i] .- layers)) + nlayers_kv[i] = k + z_layered[i] = layers[k] + end + end + else + error("""An unknown "ksat_profile" is specified in the TOML file ($ksat_profile). + This should be "exponential", "exponential_constant", "layered" or + "layered_exponential". + """) + end sbm = SBM{Float,maxlayers,maxlayers + 1}( Δt = tosecond(Δt), maxlayers = maxlayers, n = n, nlayers = nlayers, - n_unsatlayers = fill(0, n), + n_unsatlayers = n_unsatlayers, + nlayers_kv = nlayers_kv, riverfrac = riverfrac, θₛ = θₛ, θᵣ = θᵣ, kv₀ = kv₀, + kv = svectorscopy(kv, Val{maxlayers}()), kvfrac = svectorscopy(kvfrac, Val{maxlayers}()), hb = hb, soilthickness = soilthickness, act_thickl = act_thickl, - sumlayers = svectorscopy(s_layers, Val{maxlayers + 1}()), + sumlayers = sumlayers, infiltcappath = infiltcappath, infiltcapsoil = infiltcapsoil, maxleakage = maxleakage, @@ -525,6 +594,8 @@ function initialize_sbm(nc, config, riverfrac, inds) stemflow = fill(mv, n), throughfall = fill(mv, n), f = f, + z_exp = z_exp, + z_layered = z_layered, ustorelayerdepth = zero(act_thickl), satwaterdepth = satwaterdepth, zi = zi, @@ -692,6 +763,7 @@ function update_until_recharge(sbm::SBM, config) modelsnow = get(config.model, "snow", false)::Bool transfermethod = get(config.model, "transfermethod", false)::Bool ust = get(config.model, "whole_ust_available", false)::Bool # should be removed from optional setting and code? + ksat_profile = get(config.input.vertical, "ksat_profile", "exponential")::String threaded_foreach(1:sbm.n, basesize = 250) do i if modelsnow @@ -775,7 +847,7 @@ function update_until_recharge(sbm::SBM, config) # (ast) and the updated unsaturated storage for each soil layer. if transfermethod && sbm.maxlayers == 1 ustorelayerdepth = sbm.ustorelayerdepth[i][1] + infiltsoilpath - kv_z = sbm.kvfrac[i][1] * sbm.kv₀[i] * exp(-sbm.f[i] * sbm.zi[i]) + kv_z = hydraulic_conductivity_at_depth(sbm, sbm.zi[i], i, 1, ksat_profile) ustorelayerdepth, ast = unsatzone_flow_sbm( ustorelayerdepth, sbm.soilwatercapacity[i], @@ -789,7 +861,7 @@ function update_until_recharge(sbm::SBM, config) else for m = 1:n_usl l_sat = usl[m] * (sbm.θₛ[i] - sbm.θᵣ[i]) - kv_z = sbm.kvfrac[i][m] * sbm.kv₀[i] * exp(-sbm.f[i] * z[m]) + kv_z = hydraulic_conductivity_at_depth(sbm, z[m], i, m, ksat_profile) ustorelayerdepth = m == 1 ? sbm.ustorelayerdepth[i][m] + infiltsoilpath : sbm.ustorelayerdepth[i][m] + ast @@ -901,7 +973,7 @@ function update_until_recharge(sbm::SBM, config) actcapflux = 0.0 if n_usl > 0 - ksat = sbm.kvfrac[i][n_usl] * sbm.kv₀[i] * exp(-sbm.f[i] * sbm.zi[i]) + ksat = hydraulic_conductivity_at_depth(sbm, sbm.zi[i], i, n_usl, ksat_profile) ustorecapacity = sbm.soilwatercapacity[i] - satwaterdepth - sum(@view usld[1:sbm.nlayers[i]]) maxcapflux = max(0.0, min(ksat, actevapustore, ustorecapacity, satwaterdepth)) @@ -925,7 +997,13 @@ function update_until_recharge(sbm::SBM, config) actcapflux = actcapflux + toadd end end - deepksat = sbm.kvfrac[i][end] * sbm.kv₀[i] * exp(-sbm.f[i] * sbm.soilthickness[i]) + deepksat = hydraulic_conductivity_at_depth( + sbm, + sbm.soilthickness[i], + i, + sbm.nlayers[i], + ksat_profile, + ) deeptransfer = min(satwaterdepth, deepksat) actleakage = max(0.0, min(sbm.maxleakage[i], deeptransfer)) diff --git a/src/sbm_model.jl b/src/sbm_model.jl index 085fca0d4..cdac2694b 100644 --- a/src/sbm_model.jl +++ b/src/sbm_model.jl @@ -126,11 +126,15 @@ function initialize_sbm_model(config::Config) f = sbm.f .* 1000.0 zi = sbm.zi .* 0.001 soilthickness = sbm.soilthickness .* 0.001 + z_exp = sbm.z_exp .* 0.001 ssf = LateralSSF{Float}( kh₀ = kh₀, f = f, + kh = fill(mv, n), + khfrac = khfrac, zi = zi, + z_exp = z_exp, soilthickness = soilthickness, θₛ = sbm.θₛ, θᵣ = sbm.θᵣ, @@ -140,11 +144,20 @@ function initialize_sbm_model(config::Config) dw = dw, exfiltwater = fill(mv, n), recharge = fill(mv, n), - ssf = ((kh₀ .* βₗ) ./ f) .* (exp.(-f .* zi) - exp.(-f .* soilthickness)) .* dw, + ssf = fill(mv, n), ssfin = fill(mv, n), - ssfmax = ((kh₀ .* βₗ) ./ f) .* (1.0 .- exp.(-f .* soilthickness)), + ssfmax = fill(mv, n), to_river = zeros(n), ) + # update variables `ssf`, `ssfmax` and `kh` (layered profile) based on ksat_profile + ksat_profile = get(config.input.vertical, "ksat_profile", "exponential")::String + if ksat_profile == "exponential" + initialize_lateralssf_exp!(ssf::LateralSSF) + elseif ksat_profile == "exponential_constant" + initialize_lateralssf_exp_const!(ssf::LateralSSF) + elseif ksat_profile == "layered" || ksat_profile == "layered_exponential" + initialize_lateralssf_layered!(ssf::LateralSSF, sbm::SBM, ksat_profile) + end else # when the SBM model is coupled (BMI) to a groundwater model, the following # variables are expected to be exchanged from the groundwater model. @@ -368,13 +381,21 @@ end function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmModel} @unpack lateral, vertical, network, clock, config = model + ksat_profile = get(config.input.vertical, "ksat_profile", "exponential")::String + model = update_until_recharge(model) # exchange of recharge between vertical sbm concept and subsurface flow domain lateral.subsurface.recharge .= vertical.recharge ./ 1000.0 lateral.subsurface.recharge .*= lateral.subsurface.dw lateral.subsurface.zi .= vertical.zi ./ 1000.0 # update lateral subsurface flow domain (kinematic wave) - update(lateral.subsurface, network.land, network.frac_toriver) + if (ksat_profile == "layered") || (ksat_profile == "layered_exponential") + for i in eachindex(lateral.subsurface.kh) + lateral.subsurface.kh[i] = + kh_layered_profile(vertical, lateral.subsurface.khfrac[i], i, ksat_profile) + end + end + update(lateral.subsurface, network.land, network.frac_toriver, ksat_profile) model = update_after_subsurfaceflow(model) model = update_total_water_storage(model) end diff --git a/src/utils.jl b/src/utils.jl index b229be67d..11ac219c1 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -656,3 +656,173 @@ function threaded_foreach(f, x::AbstractArray; basesize::Integer) end return nothing end + +""" + hydraulic_conductivity_at_depth(sbm::SBM, z, i, ksat_profile) + +Return vertical hydraulic conductivity `kv_z` for soil layer `n` at depth `z` for vertical +concept `SBM` (at index `i`) based on hydraulic conductivity profile `ksat_profile`. +""" +function hydraulic_conductivity_at_depth(sbm::SBM, z, i, n, ksat_profile) + if ksat_profile == "exponential" + kv_z = sbm.kvfrac[i][n] * sbm.kv₀[i] * exp(-sbm.f[i] * z) + elseif ksat_profile == "exponential_constant" + if z < sbm.z_exp[i] + kv_z = sbm.kvfrac[i][n] * sbm.kv₀[i] * exp(-sbm.f[i] * z) + else + kv_z = sbm.kvfrac[i][n] * sbm.kv₀[i] * exp(-sbm.f[i] * sbm.z_exp[i]) + end + elseif ksat_profile == "layered" + kv_z = sbm.kvfrac[i][n] * sbm.kv[i][n] + elseif ksat_profile == "layered_exponential" + if z < sbm.z_layered[i] + kv_z = sbm.kvfrac[i][n] * sbm.kv[i][n] + else + n = sbm.nlayers_kv[i] + kv_z = sbm.kvfrac[i][n] * sbm.kv[i][n] * exp(-sbm.f[i] * (z - sbm.z_layered[i])) + end + else + error("unknown ksat_profile") + end + return kv_z +end + +""" + kh_layered_profile(sbm::SBM, khfrac, z, i, ksat_profile) + +Return equivalent horizontal hydraulic conductivity `kh` [m d⁻¹] for a layered soil profile +of vertical concept `SBM` (at index `i`) based on multiplication factor `khfrac` [-], water +table depth `z` [mm] and hydraulic conductivity profile `ksat_profile`. +""" +function kh_layered_profile(sbm::SBM, khfrac, i, ksat_profile) + + m = sbm.nlayers[i] + t_factor = (tosecond(basetimestep) / sbm.Δt) + if (sbm.soilthickness[i] - sbm.zi[i]) > 0.0 + transmissivity = 0.0 + sumlayers = @view sbm.sumlayers[i][2:end] + n = max(sbm.n_unsatlayers[i], 1) + + if ksat_profile == "layered" + transmissivity += (sumlayers[n] - sbm.zi[i]) * sbm.kv[i][n] + elseif ksat_profile == "layered_exponential" + if sbm.zi[i] >= sbm.z_layered[i] + zt = sbm.soilthickness[i] - sbm.z_layered[i] + j = sbm.nlayers_kv[i] + transmissivity += + sbm.kv[i][j] / sbm.f[i] * + (exp(-sbm.f[i] * (sbm.zi[i] - sbm.z_layered[i])) - exp(-sbm.f[i] * zt)) + n = m + else + transmissivity += (sumlayers[n] - sbm.zi[i]) * sbm.kv[i][n] + end + else + error("unknown ksat_profile, `layered` or `layered_exponential` are allowed") + end + n += 1 + while n <= m + if ksat_profile == "layered" + transmissivity += sbm.act_thickl[i][n] * sbm.kv[i][n] + elseif ksat_profile == "layered_exponential" + if n > sbm.nlayers_kv[i] + zt = sbm.soilthickness[i] - sbm.z_layered[i] + j = sbm.nlayers_kv[i] + transmissivity += sbm.kv[i][j] / sbm.f[i] * (1.0 - exp(-sbm.f[i] * zt)) + n = m + else + transmissivity += sbm.act_thickl[i][n] * sbm.kv[i][n] + end + end + n += 1 + end + # convert units for kh [m d⁻¹] computation (transmissivity [mm² Δt⁻¹], soilthickness + # [mm] and zi [mm]) + kh = + 0.001 * + (transmissivity / (sbm.soilthickness[i] - sbm.zi[i])) * + t_factor * + khfrac + else + if ksat_profile == "layered" + kh = 0.001 * sbm.kv[i][m] * t_factor * khfrac + elseif ksat_profile == "layered_exponential" + if sbm.zi[i] >= sbm.z_layered[i] + j = sbm.nlayers_kv[i] + kh = + 0.001 * + sbm.kv[i][j] * + exp(-sbm.f[i] * (sbm.zi[i] - sbm.z_layered[i])) * + khfrac * + t_factor + else + kh = 0.001 * sbm.kv[i][m] * t_factor * khfrac + end + else + error("unknown ksat_profile, `layered` or `layered_exponential` are allowed") + end + end + return kh +end + +"Initialize lateral subsurface variables `ssf` and `ssfmax` with `ksat_profile` `exponential`" +function initialize_lateralssf_exp!(ssf::LateralSSF) + for i in eachindex(ssf.ssf) + ssf.ssfmax[i] = + ((ssf.kh₀[i] * ssf.βₗ[i]) / ssf.f[i]) * + (1.0 - exp(-ssf.f[i] * ssf.soilthickness[i])) + ssf.ssf[i] = + ((ssf.kh₀[i] * ssf.βₗ[i]) / ssf.f[i]) * + (exp(-ssf.f[i] * ssf.zi[i]) - exp(-ssf.f[i] * ssf.soilthickness[i])) * + ssf.dw[i] + end +end + +"Initialize lateral subsurface variables `ssf` and `ssfmax` with `ksat_profile` `exponential_constant`" +function initialize_lateralssf_exp_const!(ssf::LateralSSF) + ssf_constant = @. ssf.khfrac * + ssf.kh₀ * + exp(-ssf.f * ssf.z_exp) * + ssf.βₗ * + (ssf.soilthickness - ssf.z_exp) + for i in eachindex(ssf.ssf) + ssf.ssfmax[i] = + ((ssf.khfrac[i] * ssf.kh₀[i] * ssf.βₗ[i]) / ssf.f[i]) * + (1.0 - exp(-ssf.f[i] * ssf.z_exp[i])) + ssf_constant[i] + if ssf.zi[i] < ssf.z_exp[i] + ssf.ssf[i] = + ( + ((ssf.kh₀[i] * ssf.βₗ[i]) / ssf.f[i]) * + (exp(-ssf.f[i] * ssf.zi[i]) - exp(-ssf.f[i] * ssf.z_exp[i])) + + ssf_constant[i] + ) * ssf.dw[i] + else + ssf.ssf[i] = + ssf.kh₀[i] * + exp(-ssf.f[i] * ssf.zi[i]) * + ssf.βₗ[i] * + (ssf.soilthickness[i] - ssf.zi[i]) * + ssf.dw[i] + end + end +end + +"Initialize lateral subsurface variables `ssf`, `ssfmax` and `kh` with ksat_profile` `layered` or `layered_exponential`" +function initialize_lateralssf_layered!(ssf::LateralSSF, sbm::SBM, ksat_profile) + for i in eachindex(ssf.ssf) + ssf.kh[i] = kh_layered_profile(sbm, ssf.khfrac[i], i, ksat_profile) + ssf.ssf[i] = ssf.kh[i] * (ssf.soilthickness[i] - ssf.zi[i]) * ssf.βₗ[i] * ssf.dw[i] + kh_max = 0.0 + for j = 1:sbm.nlayers[i] + if j <= sbm.nlayers_kv[i] + kh_max += sbm.kv[i][j] * sbm.act_thickl[i][j] + else + zt = sbm.soilthickness[i] - sbm.z_layered[i] + k = max(j - 1, 1) + kh_max += sbm.kv[i][k] / sbm.f[i] * (1.0 - exp(-sbm.f[i] * zt)) + break + end + end + kh_max = kh_max * ssf.khfrac[i] * 0.001 * 0.001 + ssf.ssfmax[i] = kh_max * ssf.βₗ[i] + end +end diff --git a/test/bmi.jl b/test/bmi.jl index 848393681..af6554220 100644 --- a/test/bmi.jl +++ b/test/bmi.jl @@ -20,9 +20,9 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml") @testset "model information functions" begin @test BMI.get_component_name(model) == "sbm" - @test BMI.get_input_item_count(model) == 181 - @test BMI.get_output_item_count(model) == 181 - @test BMI.get_input_var_names(model)[[1, 5, 151, 175]] == [ + @test BMI.get_input_item_count(model) == 191 + @test BMI.get_output_item_count(model) == 191 + @test BMI.get_input_var_names(model)[[1, 6, 161, 185]] == [ "vertical.nlayers", "vertical.θᵣ", "lateral.river.q", diff --git a/test/horizontal_process.jl b/test/horizontal_process.jl index ea7873af2..9e395a098 100644 --- a/test/horizontal_process.jl +++ b/test/horizontal_process.jl @@ -62,6 +62,8 @@ end 1697.05 * 1000.0, 1200.0 * 1000.0, 2443723.716252628, + 1.0, + "exponential", ), (7.410313985168225e10, 1540.1496836278836, -0.0), ), diff --git a/test/run_sbm.jl b/test/run_sbm.jl index 5f22f35a0..7e2194976 100644 --- a/test/run_sbm.jl +++ b/test/run_sbm.jl @@ -404,3 +404,109 @@ model = Wflow.run_timestep(model) @test h[501] ≈ 0.05610231297517167f0 end Wflow.close_files(model, delete_output = false) + +# test different ksat profiles +@testset "ksat profiles (SBM)" begin + i = 100 + tomlpath = joinpath(@__DIR__, "sbm_config.toml") + config = Wflow.Config(tomlpath) + config.input.vertical.kv = "kv" + config.input.vertical.z_exp = Dict("value" => 400.0) + config.input.vertical.z_layered = Dict("value" => 400.0) + + @testset "exponential profile" begin + model = Wflow.initialize_sbm_model(config) + @unpack vertical = model + z = vertical.zi[i] + kv_z = Wflow.hydraulic_conductivity_at_depth(vertical, z, i, 2, "exponential") + @test kv_z ≈ vertical.kvfrac[i][2] * vertical.kv₀[i] * exp(-vertical.f[i] * z) + @test vertical.z_exp == vertical.soilthickness + @test_throws ErrorException Wflow.kh_layered_profile( + vertical, + 100.0, + i, + "exponential", + ) + @test all(isnan.(vertical.z_layered)) + @test all(isnan.(vertical.kv[i])) + @test all(vertical.nlayers_kv .== 0) + end + + @testset "exponential constant profile" begin + config.input.vertical.ksat_profile = "exponential_constant" + model = Wflow.initialize_sbm_model(config) + @unpack vertical = model + z = vertical.zi[i] + kv_z = + Wflow.hydraulic_conductivity_at_depth(vertical, z, i, 2, "exponential_constant") + @test kv_z ≈ vertical.kvfrac[i][2] * vertical.kv₀[i] * exp(-vertical.f[i] * z) + kv_400 = Wflow.hydraulic_conductivity_at_depth( + vertical, + 400.0, + i, + 2, + "exponential_constant", + ) + kv_1000 = Wflow.hydraulic_conductivity_at_depth( + vertical, + 1000.0, + i, + 3, + "exponential_constant", + ) + @test kv_400 ≈ kv_1000 + @test_throws ErrorException Wflow.kh_layered_profile( + vertical, + 100.0, + i, + "exponential_constant", + ) + @test all(isnan.(vertical.z_layered)) + @test all(isnan.(vertical.kv[i])) + @test all(vertical.nlayers_kv .== 0) + @test all(vertical.z_exp .== 400.0) + end + + @testset "layered profile" begin + config.input.vertical.ksat_profile = "layered" + model = Wflow.initialize_sbm_model(config) + @unpack vertical = model + z = vertical.zi[i] + @test Wflow.hydraulic_conductivity_at_depth(vertical, z, i, 2, "layered") ≈ + vertical.kv[100][2] + @test Wflow.kh_layered_profile(vertical, 100.0, i, "layered") ≈ 47.508932674632355f0 + @test vertical.nlayers_kv[i] == 4 + @test vertical.z_layered == vertical.soilthickness + @test all(isnan.(vertical.z_exp)) + end + + @testset "layered exponential profile" begin + config.input.vertical.ksat_profile = "layered_exponential" + model = Wflow.initialize_sbm_model(config) + @unpack vertical = model + z = vertical.zi[i] + @test Wflow.hydraulic_conductivity_at_depth( + vertical, + z, + i, + 2, + "layered_exponential", + ) ≈ vertical.kv[i][2] + @test vertical.nlayers_kv[i] == 2 + @test Wflow.kh_layered_profile(vertical, 100.0, i, "layered_exponential") ≈ + 33.76026208801769f0 + @test all(vertical.z_layered[1:10] .== 400.0) + @test all(isnan.(vertical.z_exp)) + end + + model = Wflow.run_timestep(model) + model = Wflow.run_timestep(model) + @testset "river flow layered exponential profile" begin + q = model.lateral.river.q_av + @test sum(q) ≈ 3118.8690178033266f0 + @test q[1622] ≈ 0.000548447582354063f0 + @test q[43] ≈ 9.860543811678328f0 + end + + Wflow.close_files(model, delete_output = false) +end From faba5266553b90cecab4343cba6641853208150c Mon Sep 17 00:00:00 2001 From: Willem van Verseveld Date: Wed, 13 Mar 2024 15:14:57 +0100 Subject: [PATCH 2/2] Add floodplain discharge to inflow reservoirs and lakes (#368) * Add floodplain discarge to inflow reservoir and lake * Update changelog * Change order update reservoir and lake `ShallowWaterRiver` This way the inflow can be based on the current `q` instead of `q0` (`q` at the previous time step). Moved inflow computation of reservoir and lake to separate function for `ShallowWaterRiver`. * Fix inflow water body local inertial routing river flow For local inertial routing the inflow has been corrected: it should also include the `to_river` variable (from subsurface and overland flow), because the water body cell is not included in the local inertial routing solution (boundary condition). * Update changelog * Update Wflow ZMQ Server tests * Update tests because of merge `master` * Fix client.jl tests --- docs/src/changelog.md | 4 + server/test/client.jl | 8 +- src/flow.jl | 211 ++++++++++++++++++++----------------- test/bmi.jl | 6 +- test/horizontal_process.jl | 2 +- 5 files changed, 124 insertions(+), 107 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 318667d20..cd78cf337 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -17,6 +17,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 lake storage. The variable `actevap` has also been added to the reservoir module. - The `set_states` function for model type `sbm` with local inertial routing for river and land component. +- Inflow to reservoir and lake locations for local inertial routing: 1) with floodplain + routing, the floodplain discharge was not added to the inflow of these locations, 2) the + `to_river` variable from overland flow and lateral subsurface flow was not added to the + inflow of these locations. ### Changed - Stop exposing scalar variables through BMI. The `BMI.get_value_ptr` function was diff --git a/server/test/client.jl b/server/test/client.jl index 41fc1182d..b6e000ac8 100644 --- a/server/test/client.jl +++ b/server/test/client.jl @@ -33,15 +33,15 @@ end @testset "model information functions" begin @test request((fn = "get_component_name",)) == Dict("component_name" => "sbm") - @test request((fn = "get_input_item_count",)) == Dict("input_item_count" => 181) - @test request((fn = "get_output_item_count",)) == Dict("output_item_count" => 181) - @test request((fn = "get_input_var_names",))["input_var_names"][[1, 5, 151, 175]] == [ + @test request((fn = "get_input_item_count",)) == Dict("input_item_count" => 192) + @test request((fn = "get_output_item_count",)) == Dict("output_item_count" => 192) + @test request((fn = "get_input_var_names",))["input_var_names"][[1, 6, 161, 186]] == [ "vertical.nlayers", "vertical.θᵣ", "lateral.river.q", "lateral.river.reservoir.outflow", ] - @test request((fn = "get_output_var_names",))["output_var_names"][[1, 5, 151, 175]] == [ + @test request((fn = "get_output_var_names",))["output_var_names"][[1, 6, 161, 186]] == [ "vertical.nlayers", "vertical.θᵣ", "lateral.river.q", diff --git a/src/flow.jl b/src/flow.jl index d7dd800ac..a07f15991 100644 --- a/src/flow.jl +++ b/src/flow.jl @@ -13,6 +13,7 @@ abstract type SurfaceFlow end qlat::Vector{T} | "m2 s-1" # Lateral inflow per unit length [m² s⁻¹] inwater::Vector{T} | "m3 s-1" # Lateral inflow [m³ s⁻¹] inflow::Vector{T} | "m3 s-1" # External inflow (abstraction/supply/demand) [m³ s⁻¹] + inflow_wb::Vector{T} | "m3 s-1" # inflow waterbody (lake or reservoir model) from land part [m³ s⁻¹] volume::Vector{T} | "m3" # Kinematic wave volume [m³] (based on water level h) h::Vector{T} | "m" # Water level [m] h_av::Vector{T} | "m" # Average water level [m] @@ -158,6 +159,7 @@ function initialize_surfaceflow_river( qlat = zeros(Float, n), inwater = zeros(Float, n), inflow = zeros(Float, n), + inflow_wb = zeros(Float, n), volume = zeros(Float, n), h = zeros(Float, n), h_av = zeros(Float, n), @@ -246,7 +248,7 @@ function update(sf::SurfaceFlowLand, network, frac_toriver) sf.volume .= sf.dl .* sf.width .* sf.h end -function update(sf::SurfaceFlowRiver, network, inflow_wb, doy) +function update(sf::SurfaceFlowRiver, network, doy) @unpack graph, subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes = network @@ -306,7 +308,7 @@ function update(sf::SurfaceFlowRiver, network, inflow_wb, doy) # run reservoir model and copy reservoir outflow to inflow (qin) of # downstream river cell i = sf.reservoir_index[v] - update(sf.reservoir, i, sf.q[v] + inflow_wb[v], Δt) + update(sf.reservoir, i, sf.q[v] + sf.inflow_wb[v], Δt) downstream_nodes = outneighbors(graph, v) n_downstream = length(downstream_nodes) @@ -327,7 +329,7 @@ function update(sf::SurfaceFlowRiver, network, inflow_wb, doy) # run lake model and copy lake outflow to inflow (qin) of downstream river # cell i = sf.lake_index[v] - update(sf.lake, i, sf.q[v] + inflow_wb[v], doy, Δt) + update(sf.lake, i, sf.q[v] + sf.inflow_wb[v], doy, Δt) downstream_nodes = outneighbors(graph, v) n_downstream = length(downstream_nodes) @@ -691,14 +693,16 @@ function initialize_shallowwater_river( return sw_river, nodes_at_link end -function shallowwater_river_update( - sw::ShallowWaterRiver, - network, - Δt, - inflow_wb, - doy, - update_h, -) +"Return the upstream inflow for a waterbody in `ShallowWaterRiver`" +function get_inflow_waterbody(sw::ShallowWaterRiver, src_edge) + q_in = sum_at(sw.q, src_edge) + if !isnothing(sw.floodplain) + q_in = q_in + sum_at(sw.floodplain.q, src_edge) + end + return q_in +end + +function shallowwater_river_update(sw::ShallowWaterRiver, network, Δt, doy, update_h) @unpack nodes_at_link, links_at_node = network @@ -706,32 +710,6 @@ function shallowwater_river_update( if !isnothing(sw.floodplain) sw.floodplain.q0 .= sw.floodplain.q end - # For reservoir and lake locations the local inertial solution is replaced by the - # reservoir or lake model. These locations are handled as boundary conditions in the - # local inertial model (fixed h). - for v in eachindex(sw.reservoir_index) - i = sw.reservoir_index[v] - update( - sw.reservoir, - v, - sum_at(sw.q0, links_at_node.src[i]) + inflow_wb[i] + sw.inflow_wb[i], - Δt, - ) - sw.q[i] = sw.reservoir.outflow[v] - sw.q_av[i] += sw.q[i] * Δt - end - for v in eachindex(sw.lake_index) - i = sw.lake_index[v] - update( - sw.lake, - v, - sum_at(sw.q0, links_at_node.src[i]) + inflow_wb[i] + sw.inflow_wb[i], - doy, - Δt, - ) - sw.q[i] = sw.lake.outflow[v] - sw.q_av[i] += sw.q[i] * Δt - end @tturbo for j in eachindex(sw.active_e) i = sw.active_e[j] i_src = nodes_at_link.src[i] @@ -861,6 +839,25 @@ function shallowwater_river_update( sw.floodplain.q_av[i] += sw.floodplain.q[i] * Δt end end + # For reservoir and lake locations the local inertial solution is replaced by the + # reservoir or lake model. These locations are handled as boundary conditions in the + # local inertial model (fixed h). + for v in eachindex(sw.reservoir_index) + i = sw.reservoir_index[v] + + q_in = get_inflow_waterbody(sw, links_at_node.src[i]) + update(sw.reservoir, v, q_in + sw.inflow_wb[i], Δt) + sw.q[i] = sw.reservoir.outflow[v] + sw.q_av[i] += sw.q[i] * Δt + end + for v in eachindex(sw.lake_index) + i = sw.lake_index[v] + + q_in = get_inflow_waterbody(sw, links_at_node.src[i]) + update(sw.lake, v, q_in + sw.inflow_wb[i], doy, Δt) + sw.q[i] = sw.lake.outflow[v] + sw.q_av[i] += sw.q[i] * Δt + end if update_h @batch per = thread minbatch = 2000 for i in sw.active_n @@ -908,13 +905,7 @@ function shallowwater_river_update( end end -function update( - sw::ShallowWaterRiver{T}, - network, - inflow_wb, - doy; - update_h = true, -) where {T} +function update(sw::ShallowWaterRiver{T}, network, doy; update_h = true) where {T} @unpack nodes_at_link, links_at_node = network if !isnothing(sw.reservoir) @@ -940,7 +931,7 @@ function update( if t + Δt > sw.Δt Δt = sw.Δt - t end - shallowwater_river_update(sw, network, Δt, inflow_wb, doy, update_h) + shallowwater_river_update(sw, network, Δt, doy, update_h) t = t + Δt end sw.q_av ./= sw.Δt @@ -990,6 +981,7 @@ const dirs = (:yd, :xd, :xu, :yu) volume::Vector{T} | "m3" # total volume of cell (including river volume for river cells) error::Vector{T} | "m3" # error volume runoff::Vector{T} | "m3 s-1" # runoff from hydrological model + inflow_wb::Vector{T} | "m3 s-1" # inflow to water body from hydrological model h::Vector{T} | "m" # water depth of cell (for river cells the reference is the river bed elevation `zb`) z::Vector{T} | "m" # elevation of cell froude_limit::Bool | "-" | 0 | "none" | "none" # if true a check is performed if froude number > 1.0 (algorithm is modified) @@ -1106,6 +1098,7 @@ function initialize_shallowwater_land( volume = zeros(n), error = zeros(n), runoff = zeros(n), + inflow_wb = zeros(n), h = zeros(n), h_av = zeros(n), z = elevation, @@ -1152,7 +1145,6 @@ function update( sw::ShallowWaterLand{T}, swr::ShallowWaterRiver{T}, network, - inflow_wb, doy; update_h = false, ) where {T} @@ -1181,8 +1173,8 @@ function update( if t + Δt > swr.Δt Δt = swr.Δt - t end - shallowwater_river_update(swr, network.river, Δt, inflow_wb, doy, update_h) - update(sw, swr, network, Δt) + shallowwater_river_update(swr, network.river, Δt, doy, update_h) + shallowwater_update(sw, swr, network, Δt) t = t + Δt end swr.q_av ./= swr.Δt @@ -1190,7 +1182,12 @@ function update( sw.h_av ./= sw.Δt end -function update(sw::ShallowWaterLand{T}, swr::ShallowWaterRiver{T}, network, Δt) where {T} +function shallowwater_update( + sw::ShallowWaterLand{T}, + swr::ShallowWaterRiver{T}, + network, + Δt, +) where {T} indices = network.land.staggered_indices inds_riv = network.land.index_river @@ -1296,7 +1293,9 @@ function update(sw::ShallowWaterLand{T}, swr::ShallowWaterRiver{T}, network, Δt # for reservoir or lake set inflow from land part, these are boundary points # and update of volume and h is not required swr.inflow_wb[inds_riv[i]] = - sw.runoff[i] + (sw.qx[xd] - sw.qx[i] + sw.qy[yd] - sw.qy[i]) + sw.inflow_wb[i] + + sw.runoff[i] + + (sw.qx[xd] - sw.qx[i] + sw.qy[yd] - sw.qy[i]) else sw.volume[i] += ( @@ -1587,7 +1586,7 @@ end """ set_river_inwater(model::Model{N,L,V,R,W,T}, ssf_toriver) where {N,L,V<:SBM,R,W,T} -Set `inwater` of the river component for a `Model` with vertical `SBM` concept. +Set `inwater` of the lateral river component for a `Model` with vertical `SBM` concept. `ssf_toriver` is the subsurface flow to the river. """ function set_river_inwater(model::Model{N,L,V,R,W,T}, ssf_toriver) where {N,L,V<:SBM,R,W,T} @@ -1612,7 +1611,7 @@ end """ set_river_inwater(model, ssf_toriver) -Set `inwater` of the river component (based on overland flow). +Set `inwater` of the lateral river component (based on overland flow). """ function set_river_inwater(model, ssf_toriver) @unpack lateral, network = model @@ -1623,7 +1622,7 @@ end """ set_land_inwater(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel} -Set `inwater` of the land component for the `SbmGwgModel` type. +Set `inwater` of the lateral land component for the `SbmGwgModel` type. """ function set_land_inwater(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel} @unpack lateral, vertical, network, config = model @@ -1643,7 +1642,7 @@ end """ set_land_inwater(model) -Set `inwater` of the land component, based on `runoff` of the `vertical` concept. +Set `inwater` of the lateral land component, based on `runoff` of the `vertical` concept. """ function set_land_inwater(model) @unpack lateral, vertical, network = model @@ -1651,66 +1650,90 @@ function set_land_inwater(model) (vertical.runoff .* network.land.xl .* network.land.yl .* 0.001) ./ lateral.land.Δt end +# Computation of inflow from the lateral components `land` and `subsurface` to water bodies +# depends on the routing scheme (see different `get_inflow_waterbody` below). For the river +# kinematic wave, the variables `to_river` can be excluded, because this part is added to +# the river kinematic wave (kinematic wave is also solved for the water body cell). For +# local inertial river routing, `to_river` is included, because for the local inertial +# solution the water body cells are excluded (boundary condition). For `GroundwaterFlow` +# (Darcian flow in 4 directions), the lateral subsurface flow is excluded (for now) and +# inflow consists of overland flow. """ - get_inflow_waterbody(model) + set_inflow_waterbody( + model::Model{N,L,V,R,W,T}, + ) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,SurfaceFlow,SurfaceFlow}},V,R,W,T} -Get inflow to a water body (reservoir or lake) `inflow_wb` based on overland flow. +Set inflow from the subsurface and land components to a water body (reservoir or lake) +`inflow_wb` from a model type that contains the lateral components `SurfaceFlow`. """ -function get_inflow_waterbody(model) +function set_inflow_waterbody( + model::Model{N,L,V,R,W,T}, +) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,SurfaceFlow,SurfaceFlow}},V,R,W,T} @unpack lateral, network = model + @unpack subsurface, land, river = lateral inds = network.index_river + if !isnothing(lateral.river.reservoir) || !isnothing(lateral.river.lake) - inflow_wb = lateral.land.q_av[inds] - else - inflow_wb = nothing + if typeof(subsurface) <: LateralSSF || typeof(subsurface) <: GroundwaterExchange + @. river.inflow_wb = + subsurface.ssf[inds] / tosecond(basetimestep) + land.q_av[inds] + elseif typof(subsurface.flow) <: GroundwaterFlow || isnothing(subsurface) + river.inflow_wb .= land.q_av[inds] + end end - return inflow_wb end """ - get_inflow_waterbody( + set_inflow_waterbody( model::Model{N,L,V,R,W,T}, - ) where {N,L<:NamedTuple{<:Any,<:Tuple{LateralSSF,SurfaceFlow,Any}},V,R,W,T} + ) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,SurfaceFlow,ShallowWaterRiver}},V,R,W,T} -Get inflow to a water body (reservoir or lake) `inflow_wb` from a model type that contains -the lateral components `LateralSSF` and `SurfaceFlow`. +Set inflow from the subsurface and land components to a water body (reservoir or lake) +`inflow_wb` from a model type that contains the lateral components `SurfaceFlow` and +`ShallowWaterRiver`. """ -function get_inflow_waterbody( +function set_inflow_waterbody( model::Model{N,L,V,R,W,T}, -) where {N,L<:NamedTuple{<:Any,<:Tuple{LateralSSF,SurfaceFlow,Any}},V,R,W,T} +) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,SurfaceFlow,ShallowWaterRiver}},V,R,W,T} @unpack lateral, network = model - + @unpack subsurface, land, river = lateral inds = network.index_river + if !isnothing(lateral.river.reservoir) || !isnothing(lateral.river.lake) - inflow_wb = - lateral.subsurface.ssf[inds] ./ tosecond(basetimestep) .+ - lateral.land.q_av[inds] - else - inflow_wb = nothing + if typeof(subsurface) <: LateralSSF || typeof(subsurface) <: GroundwaterExchange + @. river.inflow_wb = + (subsurface.ssf[inds] + subsurface.to_river[inds]) / + tosecond(basetimestep) + + land.q_av[inds] + + land.to_river[inds] + elseif typeof(subsurface.flow) <: GroundwaterFlow || isnothing(subsurface) + @. river.inflow_wb = lateral.land.q_av[inds] + lateral.land.to_river[inds] + end end - return inflow_wb end """ - get_inflow_waterbody( + set_inflow_waterbody( model::Model{N,L,V,R,W,T}, - ) where {N,L<:NamedTuple{<:Any,<:Tuple{LateralSSF,ShallowWaterLand,Any}},V,R,W,T} + ) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,ShallowWaterLand,ShallowWaterRiver}},V,R,W,T} -Get inflow to a water body (reservoir or lake) `inflow_wb` from a model type that contains -the lateral components `LateralSSF` and `ShallowWaterLand`. +Set inflow from the subsurface and land components to a water body (reservoir or lake) +`inflow_wb` from a model type that contains the lateral components `ShallowWaterLand` and +`ShallowWaterRiver`. """ -function get_inflow_waterbody( +function set_inflow_waterbody( model::Model{N,L,V,R,W,T}, -) where {N,L<:NamedTuple{<:Any,<:Tuple{LateralSSF,ShallowWaterLand,Any}},V,R,W,T} +) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,ShallowWaterLand,ShallowWaterRiver}},V,R,W,T} @unpack lateral, network = model - + @unpack subsurface, land, river = lateral inds = network.index_river + if !isnothing(lateral.river.reservoir) || !isnothing(lateral.river.lake) - inflow_wb = lateral.subsurface.ssf[inds] ./ tosecond(basetimestep) - else - inflow_wb = nothing + if typeof(subsurface) <: LateralSSF || typeof(subsurface) <: GroundwaterExchange + @. land.inflow_wb[inds] = + (subsurface.ssf[inds] + subsurface.to_river[inds]) / tosecond(basetimestep) + end end - return inflow_wb end """ @@ -1728,12 +1751,8 @@ function surface_routing(model; ssf_toriver = 0.0) # run river flow set_river_inwater(model, ssf_toriver) - update( - lateral.river, - network.river, - get_inflow_waterbody(model), - julian_day(clock.time - clock.Δt), - ) + set_inflow_waterbody(model) + update(lateral.river, network.river, julian_day(clock.time - clock.Δt)) end """ @@ -1761,12 +1780,6 @@ function surface_routing( vertical.Δt ) ) - - update( - lateral.land, - lateral.river, - network, - get_inflow_waterbody(model), - julian_day(clock.time - clock.Δt), - ) + set_inflow_waterbody(model) + update(lateral.land, lateral.river, network, julian_day(clock.time - clock.Δt)) end diff --git a/test/bmi.jl b/test/bmi.jl index af6554220..e972574fb 100644 --- a/test/bmi.jl +++ b/test/bmi.jl @@ -20,9 +20,9 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml") @testset "model information functions" begin @test BMI.get_component_name(model) == "sbm" - @test BMI.get_input_item_count(model) == 191 - @test BMI.get_output_item_count(model) == 191 - @test BMI.get_input_var_names(model)[[1, 6, 161, 185]] == [ + @test BMI.get_input_item_count(model) == 192 + @test BMI.get_output_item_count(model) == 192 + @test BMI.get_input_var_names(model)[[1, 6, 161, 186]] == [ "vertical.nlayers", "vertical.θᵣ", "lateral.river.q", diff --git a/test/horizontal_process.jl b/test/horizontal_process.jl index 9e395a098..851429dd7 100644 --- a/test/horizontal_process.jl +++ b/test/horizontal_process.jl @@ -232,7 +232,7 @@ end sw_river.inwater[1] = 20.0 h0 = mean(sw_river.h) Δt = Wflow.stable_timestep(sw_river) - Wflow.shallowwater_river_update(sw_river, network, Δt, nothing, 0.0, true) + Wflow.shallowwater_river_update(sw_river, network, Δt, 0.0, true) d = abs(h0 - mean(sw_river.h)) if d <= ϵ break