{ "cells": [ { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "# 2. Slug Test - Falling Head\n", "**This test is taken from examples of AQTESOLV.**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 1. Import required libraries" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "\n", "import ttim as ttm\n", "\n", "plt.rcParams[\"figure.figsize\"] = [5, 3]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction and Conceptual Model\n", "\n", "In this notebook, we reproduce the work of Yang (2020) to check the TTim performance in analysing slug-test. We later compare the solution in TTim with the KGS analytical model (Hyder et al. 1994) implemented in AQTESOLV (Duffield, 2007).\n", "\n", "This slug test was reported in Batu (1998). A well partially penetrates a sandy unconfined aquifer that has a saturated depth of 32.57 ft. The top of the screen is located 0.47 ft below the water table and has 13.8 ft in length. The well and casing radii are 5 and 2 inches, respectively.\n", "\n", "The slug displacement is 1.48 ft. Head change has been recorded at the slug well.\n", "\n", "The conceptual model is seen in the figure below." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "editable": true, "jupyter": { "source_hidden": true }, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdsAAAE6CAYAAABeYvgTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJvklEQVR4nO3dd1hT598G8DsgJGxBtspyKyqKVkGr4sBJHT9XrRbUWvce1dYquHBr1bpna121DtTairvWUbRYV90oqCCKCogKSJ73D19SYwImmBDG/bmuXJrnPDm5T07IN+fkOedIhBACREREpDdGhg5ARERU1LHYEhER6RmLLRERkZ6x2BIREekZiy0REZGesdgSERHpGYstERGRnrHYEhER6RmLLRERkZ6x2OazCxcuoHfv3vD09IRMJoOlpSVq166N2bNn48mTJ4aOp3MPHjxAaGgozp8/ny/P16RJEzRp0kSjfhKJBF5eXlB3ErXjx49DIpFAIpFg/fr1Osu3fv16SCQS3LlzR+vHhoaGQiKR6CwLANy5c0exnO/e6tSpo9W8PDw8EBISojLvt1+/D1n+D5X93GfPnlU7vV27dvDw8MjfUP9P03UbEhKS4/rS9XvDUNS9b4qCEoYOUJysWrUKgwYNQqVKlTB27FhUrVoVmZmZOHv2LJYvX45Tp05h586dho6pUw8ePEBYWBg8PDzg4+Nj6DhKrKysEBMTg8OHD6NZs2ZK09auXQtra2ukpKQYKF3+Gjp0KHr06KHUZmlpqfPnadu2LU6dOgUXFxedz7u4MDMzw+HDhw0dg7TEYptPTp06hYEDB6JFixbYtWsXpFKpYlqLFi0wevRo/PbbbwZMWPy4ubnBysoKa9euVSq2qamp+Pnnn/HZZ59h1apVBkyYf9zc3FC/fn29P4+DgwMcHBz0/jxFmZGRUb6sK9It7kbOJzNmzIBEIsHKlSuVCm02U1NTfPLJJ4r7crkcs2fPRuXKlSGVSuHo6IjPP/8c9+7dU3pckyZN4O3tjaioKHz88ccwNzeHl5cXZs6cCblcrtT32bNnGD16NLy8vBTzbNOmDa5evarok5GRgWnTpime18HBAb1798ajR4+U5uXh4YF27dph586dqFGjBmQyGby8vLBo0SJFn6NHj6Ju3boAgN69eyt2dYWGhiqyq9vlGxISorI7LywsDPXq1YOdnR2sra1Ru3ZtrFmzRu0uYG306dMHO3bswLNnzxRtW7ZsAQB0795d7WNOnDiBZs2awcrKCubm5vD398e+fftU+p0+fRoNGjSATCaDq6srJkyYgMzMTLXz3Lp1K/z8/GBhYQFLS0u0bNkS0dHRH7RsuvDq1SuMHj0aPj4+sLGxgZ2dHfz8/LB79+48zU/dbmRt3sOXL19GYGAgzM3N4eDggMGDB2Pfvn2QSCQ4evToByypekIILF26FD4+PjAzM4OtrS06d+6M27dvK/WLjIxE+/btUaZMGchkMpQvXx79+/fH48ePVea5b98++Pj4QCqVwtPTE3PnztV57gEDBkAmk+HcuXOKNrlcjmbNmsHJyQnx8fEAgEePHmHQoEGoWrUqLC0t4ejoiKZNm+KPP/5Qml/2rt05c+Zg1qxZ8PDwgJmZGZo0aYLr168jMzMT48ePh6urK2xsbNCxY0ckJiYqzUOTz4zc3LhxAz169ICjoyOkUimqVKmC77///gNfqfzDYpsPsrKycPjwYfj6+qJs2bIaPWbgwIH46quv0KJFC0RERGDq1Kn47bff4O/vr/IHnJCQgM8++ww9e/ZEREQEWrdujQkTJmDjxo2KPqmpqWjYsCFWrFiB3r17Y8+ePVi+fDkqVqyo+MOTy+Vo3749Zs6ciR49emDfvn2YOXMmIiMj0aRJE7x8+VLpec+fP48RI0Zg5MiR2LlzJ/z9/TF8+HDFh0ft2rWxbt06AMDEiRNx6tQpnDp1Cl988YXWr+GdO3fQv39/bNu2DTt27ECnTp0wdOhQTJ06Vet5va179+4wNjbG5s2bFW1r1qxB586dYW1trdL/2LFjaNq0KZKTk7FmzRps3rwZVlZWCAoKwtatWxX9rly5gmbNmuHZs2dYv349li9fjujoaEybNk1lnjNmzMCnn36KqlWrYtu2bfjxxx+RmpqKjz/+GFeuXPmg5dOUXC7H69evlW5CCKSnp+PJkycYM2YMdu3ahc2bN6Nhw4bo1KkTfvjhB509vybv4fj4eDRu3BjXrl3DsmXL8MMPPyA1NRVDhgzR6rmysrJUljV7ed/Vv39/jBgxAs2bN8euXbuwdOlSXL58Gf7+/nj48KGi361bt+Dn54dly5bhwIEDmDRpEs6cOYOGDRsqfcE6dOgQ2rdvDysrK2zZsgVz5szBtm3bFH8nmlKX/+0vJgsXLkSVKlXQtWtXxRfJsLAwHD16FBs3blTsxs8eJzJ58mTs27cP69atg5eXF5o0aaL2y8v333+PP//8E99//z1Wr16Nq1evIigoCH379sWjR4+wdu1azJ49GwcPHlT7d/6+z4ycXLlyBXXr1sWlS5cwb9487N27F23btsWwYcMQFham1WtnMIL0LiEhQQAQ3bt316j/v//+KwCIQYMGKbWfOXNGABBff/21oq1x48YCgDhz5oxS36pVq4qWLVsq7k+ZMkUAEJGRkTk+7+bNmwUA8csvvyi1R0VFCQBi6dKlijZ3d3chkUjE+fPnlfq2aNFCWFtbi7S0NKXHrlu3TuX5GjduLBo3bqzSHhwcLNzd3XPMmZWVJTIzM8WUKVNEqVKlhFwuf+881T13tWrVFM9Xp04dIYQQly9fFgDE0aNH1WavX7++cHR0FKmpqYq2169fC29vb1GmTBlFlm7dugkzMzORkJCg1K9y5coCgIiJiRFCCBEbGytKlCghhg4dqpQvNTVVODs7i65duyraJk+eLHT9JxsTEyMAqL2pe6+8fv1aZGZmir59+4patWopTXN3dxfBwcEq83779Vu3bp3S8guh+Xt47NixQiKRiMuXLyv1a9mypQAgjhw5kuuyZj93bre333enTp0SAMS8efOU5hMXFyfMzMzEuHHj1D6PXC4XmZmZ4u7duwKA2L17t2JavXr1hKurq3j58qWiLSUlRdjZ2Wm0boODg3PM3qxZM6W+N27cENbW1qJDhw7i4MGDwsjISEycODHX+Wev32bNmomOHTsq2rPXZc2aNUVWVpaifeHChQKA+OSTT5TmM2LECAFAJCcnK9o0/cxQ975p2bKlKFOmjNL8hBBiyJAhQiaTiSdPnuS6XAUBt2wLoCNHjgCA0shOAPjoo49QpUoVHDp0SKnd2dkZH330kVJbjRo1cPfuXcX9/fv3o2LFimjevHmOz7t3716ULFkSQUFBSt+YfXx84OzsrPJNt1q1aqhZs6ZSW48ePZCSkoK///5b08XVyOHDh9G8eXPY2NjA2NgYJiYmmDRpEpKSklR2V2mrT58+OHv2LC5evIg1a9agXLlyaNSokUq/tLQ0nDlzBp07d1YaPGRsbIxevXrh3r17uHbtGoA36zB7l93b/bp166Y0z99//x2vX7/G559/rvSay2QyNG7cWOtdo0IIlS0eTQwfPhxRUVFKt3r16gEAfv75ZzRo0ACWlpYoUaIETExMsGbNGvz7779aZcuNJu/hY8eOwdvbG1WrVlXq9+mnn2r1XD/88IPKskZFRaFhw4ZK/fbu3QuJRIKePXsqvZ7Ozs6oWbOm0rpJTEzEgAEDULZsWcVr5O7uDgCK1yktLQ1RUVHo1KkTZDKZ4rHZe0Y0ZWZmpjb/0qVLlfqVL18eq1atwq5du9CuXTt8/PHHip9w3rZ8+XLUrl0bMplMkf3QoUNq12+bNm1gZPRf2ahSpQqANwPf3pbdHhsbq9Sel8+MV69e4dChQ+jYsSPMzc2V1kWbNm3w6tUrnD59OodXq+DgAKl8YG9vD3Nzc8TExGjUPykpCQDUjth0dXVV+gACgFKlSqn0k0qlSrt9Hz16BDc3t1yf9+HDh3j27BlMTU3VTn9397Wzs7NKn+y27GXQhb/++guBgYFo0qQJVq1ahTJlysDU1BS7du3C9OnTVXZva6tRo0aoUKECVqxYgW3btmHEiBFqD6N4+vQphBA5rhfgv+VOSkrK9fXJlr0rMvu37Xe9/cGmiWPHjiEgIECpLSYm5r2HtJQpU0btoT47duxA165d0aVLF4wdOxbOzs4oUaIEli1bhrVr12qVLTeavIeTkpLg6emp0u/tLzSaqFKlitpltbGxQVxcnOL+w4cPIYTIcf5eXl4A3uyCDwwMxIMHD/Dtt9+ievXqsLCwgFwuR/369RXL8PTpU8jlco3eF7kxMjLS+LCstm3bwsnJCQ8fPsSoUaNgbGysNH3+/PkYPXo0BgwYgKlTp8Le3h7Gxsb49ttv1RZbOzs7pfvZnxU5tb969UqpPS+fGUlJSXj9+jUWL16MxYsXq+2j7rfxgobFNh8YGxujWbNm2L9/P+7du4cyZcrk2j/7gyc+Pl6l74MHD2Bvb691BgcHB5XBVe+yt7dHqVKlchwVbWVlpXQ/ISFBpU92m7oPz3fJZDIkJyertL/7h7NlyxaYmJhg7969SlsEu3bteu9zaKp3796YOHEiJBIJgoOD1faxtbWFkZGR4jfutz148AAAFOumVKlSub4+2bL7b9++XbEl9CF8fX0RFRWl1Jb9RSAvNm7cCE9PT2zdulXpC0h6enqe55lXpUqVUvqdNJu611kX7O3tIZFI8Mcff6gd1JjddunSJfzzzz9Yv3690nvn5s2bSv1tbW0hkUg0el/oyoABA5Camopq1aph2LBh+Pjjj2Fra6uYvnHjRjRp0gTLli1Telxqaqpe8uTlM8PW1lax92jw4MFq+6j7ElbQcDdyPpkwYQKEEOjXrx8yMjJUpmdmZmLPnj0AgKZNmwKA0uAQAIiKisK///6rckyoJlq3bo3r16/nenxeu3btkJSUhKysLNSpU0flVqlSJaX+ly9fxj///KPUtmnTJlhZWaF27doA/vtAUrf16eHhgevXryt9cCclJeHkyZNK/SQSCUqUKKH0rfzly5f48ccfNVz69wsODkZQUBDGjh2L0qVLq+1jYWGBevXqYceOHUrLI5fLsXHjRpQpUwYVK1YEAAQEBODQoUNKxSErK0tpEBUAtGzZEiVKlMCtW7fUvubanljCyspK5fE57anQhEQigampqVKhTUhIyPNo5A/RuHFjXLp0SWXQWPbocV1r164dhBC4f/++2vVSvXp1AFC8Nu8W5BUrVijdt7CwwEcffYQdO3YobfGlpqYq/vZ1afXq1di4cSOWLFmCiIgIPHv2DL1791bqI5FIVHJfuHABp06d0nkeQLPPjHeZm5sjICAA0dHRqFGjhtp1ocmXe0Pjlm0+yR6pOGjQIPj6+mLgwIGoVq0aMjMzER0djZUrV8Lb2xtBQUGoVKkSvvzySyxevBhGRkZo3bo17ty5g2+//RZly5bFyJEjtX7+ESNGYOvWrWjfvj3Gjx+Pjz76CC9fvsSxY8fQrl07BAQEoHv37vjpp5/Qpk0bDB8+HB999BFMTExw7949HDlyBO3bt0fHjh0V83R1dcUnn3yC0NBQuLi4YOPGjYiMjMSsWbNgbm4OAChXrhzMzMzw008/oUqVKrC0tISrqytcXV3Rq1cvrFixAj179kS/fv2QlJSE2bNnq4wCbtu2LebPn48ePXrgyy+/RFJSEubOnat2ayOvXF1dNdpSDg8PR4sWLRAQEIAxY8bA1NQUS5cuxaVLl7B582bFB+/EiRMRERGBpk2bYtKkSTA3N8f333+PtLQ0pfl5eHhgypQp+Oabb3D79m20atUKtra2ePjwIf766y9YWFgYdLRlu3btsGPHDgwaNAidO3dGXFwcpk6dChcXF9y4cSNfs4wYMQJr165F69atMWXKFDg5OWHTpk2KQ9e03eX+Pg0aNMCXX36J3r174+zZs2jUqBEsLCwQHx+PEydOoHr16hg4cCAqV66McuXKYfz48RBCwM7ODnv27EFkZKTKPKdOnYpWrVopjq3PysrCrFmzYGFhofEZ5ORyeY6/UdaqVQtSqRQXL17EsGHDEBwcrCiw2aPsFy5ciBEjRgB4s36nTp2KyZMnK0Z6T5kyBZ6enhr/3q8NTT4z1Pnuu+/QsGFDfPzxxxg4cCA8PDyQmpqKmzdvYs+ePYXjJB+GHJ1VHJ0/f14EBwcLNzc3YWpqKiwsLEStWrXEpEmTRGJioqJfVlaWmDVrlqhYsaIwMTER9vb2omfPniIuLk5pfm+Pqn2buhG9T58+FcOHDxdubm7CxMREODo6irZt24qrV68q+mRmZoq5c+eKmjVrCplMJiwtLUXlypVF//79xY0bNxT93N3dRdu2bcX27dtFtWrVhKmpqfDw8BDz589XybJ582ZRuXJlYWJiIgCIyZMnK6Zt2LBBVKlSRchkMlG1alWxdetWtdnXrl0rKlWqJKRSqfDy8hLh4eFizZo1ake2ajsaOSc5jaT+448/RNOmTYWFhYUwMzMT9evXF3v27FF5/J9//inq168vpFKpcHZ2FmPHjhUrV65UySyEELt27RIBAQHC2tpaSKVS4e7uLjp37iwOHjyo6KPP0chz5szJsc/MmTOFh4eHkEqlokqVKmLVqlVqs3zIaGRN38OXLl0SzZs3FzKZTNjZ2Ym+ffuKDRs2CADin3/+yXVZs587KipK7fS2bduqHQW/du1aUa9ePcX6LleunPj888/F2bNnFX2uXLkiWrRoIaysrIStra3o0qWLiI2NVXm/CyFERESEqFGjhjA1NRVubm5i5syZGq/b3EYjAxA3btwQz58/F5UrVxZVq1ZVjPDNNnjwYGFiYqIY+Z2eni7GjBkjSpcuLWQymahdu7bYtWuXymuf0/vkyJEjAoD4+eefldrVvdaafmaoe99kt/fp00eULl1amJiYCAcHB+Hv7y+mTZv23tetIJAI8YFnBaBiycPDA97e3ti7d6+ho1Ax9+WXX2Lz5s1ISkr6oF3mpF/F/TODu5GJqNCYMmUKXF1d4eXlhefPn2Pv3r1YvXo1Jk6cyEJLBRqLLREVGiYmJpgzZw7u3buH169fo0KFCpg/fz6GDx9u6GhEueJuZCIiIj0rNIf+hIeHo27durCysoKjoyM6dOigOFtPNiEEQkND4erqqjhJ9uXLlw2UmIiI6I1CU2yPHTuGwYMH4/Tp04iMjMTr168RGBiodCjF7NmzMX/+fCxZsgRRUVFwdnZGixYt9HaANhERkSYK7W7kR48ewdHREceOHUOjRo0ghICrqytGjBiBr776CsCbs9w4OTlh1qxZ6N+/v4ETExFRcVVoB0hln+Yv+5ycMTExSEhIQGBgoKKPVCpF48aNcfLkyRyLbXp6utIZjORyOZ48eYJSpUqpPT8uEREVfUIIpKamwtXVVScnTCmUxVYIgVGjRqFhw4bw9vYG8N/5Nd89abiTk5PKifvfFh4eXniuh0hERPkqLi7uveez10ShLLZDhgzBhQsXcOLECZVp726NCiFy3UKdMGECRo0apbifnJwMNzc3nD0bB0tL1YuHExFR0ff8eQrq1CmrcgGWvCp0xXbo0KGIiIjA8ePHlb5tZF+mKSEhQekSaImJiblegksqlao9x66lpTWsrFhsiYiKM139nFhoRiMLITBkyBDs2LEDhw8fVrmkkqenJ5ydnZVO/p2RkYFjx47B398/v+MSEREpFJot28GDB2PTpk3YvXs3rKysFL/R2tjYwMzMDBKJBCNGjMCMGTNQoUIFVKhQATNmzIC5uTl69Ohh4PRERFScFZpim31x4yZNmii1r1u3DiEhIQCAcePG4eXLlxg0aBCePn2KevXq4cCBAzrb505ERJQXhfY4W31JSUmBjY0Nrl5N5m+2RETFVGpqCipXtkFycrLKNbbzotD8ZktERFRYsdgSERHpGYstERGRnrHYEhER6RmLLRERkZ6x2BIREekZiy0REZGesdgSERHpGYstERGRnrHYEhER6RmLLREp2fbzcfg3HIV//rlt6ChERQaLLREpvHjxClOnb8bd2ERMmbYJPHU6kW6w2BKRwsrVv+Hp0+cAgNNnruLI0Qtq+x069CtKl5Yobm5uJVCtWil06tQIq1d/h/T0dL3mfP48FdOmjcOnnwaienUHlC4twbx5oXme36ZNq1G6tAQVKliqTLt0KRp9+nRA7dquKFfOHI0aVcaCBVPw8uULjea9e/dWBARUQ7lyZihdWoJLl84jKuok5s0LRXLyszxnpsKFxZaIAACPHydj8ZIIxdaskZEEYVN+QlaWXKXvxYt/AwBWr96BiIhT2L79KGbPXony5atg6tQxCAqqr9dC8vRpEn76aSUyMtLRqlWHD5pXfPx9TJ06Bs7OrirTrl+/gvbt/REXdwehoQuxYcNetG/fHQsWTMGgQZ++d95JSY8wfHgvuLuXw8aNvyEi4hTKlauIc+dOYv78MKSkPPug7FR4FJrr2RKRfi34bhcyM18r7svlAjdvPcD2X/5At66NlfpeuvQ3LCws0apVB0gkEkV727b/Q506/hg5MgRz507C1KmL9JK1TBl3XLnyFBKJBE+ePMamTavzPK/x4wegXr1GKFnSDvv2bVeatnPnJrx69QqrVv0CD49yAICGDZvi4cN4/PTTSjx79hQlS9rmOO/bt68jMzMTnTr1hJ9f4xz7UdHHLVsiQkbGa2zZegxyuUCJEsYAoPh3/YaDKv0vXDiHatV8lApttq5dg+HiUga//vqL3vJKJBK1z62tX37ZiNOnj2HGjKVqp5uYmAAArK1tlNptbErCyMgIpqamOc57xIgQdOjQEAAwcGA3lC4tQefOTTBvXiimTh0LAKhf31OxK/7kyaMfvDxUcHHLlohgaloC8+d+iciDf2PnrpMAgNevs/Blv9ZoFeir1PfJkyTcvx+Lli3b5zg/L6+KOHnyCORyOYyMlL/TCyGQlZWlUa4SJfT3EfX4cSImTx6BCRNmwtW1jNo+XboEY/XqhRg/fiC++WYWSpVywKlTx7Bx4wqEhAyGublFjvMfMeJb+Ph8hG++GYzx42fA3z8AVlbWsLS0xrNnT7B27WKsXr0Djo4uAICKFavqZTmpYGCxJSIAQPtP6sPOzlJRbAFg+NAOKFlSuaBcuvTm91pv71o5zisjIx1mZuYqhRYATp06hi5dAjTKdPp0DMqW9dCor7YmTBiEcuUqITh4YI59ypb1QETEKfTt2xH+/uUU7X37DkNY2MJc5+/hUU5RQD09K8DXt75iWunSbgDevIb6Wj4qWFhsiUgrFy6cAwBUq+aTY587d26iXLlKaqfVqOGLX3+N0ui5nJxUBy3pwr59v+DgwT34/ffoXHdHx8XdQXBwEBwcnLBy5XaUKuWA6Ogz+O67aUhLe45589boJR8VPSy2RKSVixf/homJCSpWrKZ2+vnzUXj06CE+++xLtdMtLCxzLdRv08du5LS05/jmm8Ho3XsonJxcFaOmMzMzAADJyc9gYmICc3MLzJgxHs+fpyAy8rxil3H9+o1gZ2ePUaP6oHPnzznwiTTCYktEWrl06W9UrFhN7eAgIQRmz54ImcwMwcGD1D7e0LuRnzx5jEePHmLFinlYsWKeyvSqVW3RsmV7rF27C5cvn0eFClVVfputWbMuAODatUsstqQRFlsi0lhKSjLu3r2Nrl1DVKZlZmZi4sShOHbsAGbOXA5HR2e18zD0bmQHB2f8/PMRlfbvv5+J06eP4ccf98POzl7x/NeuXUJa2nNYWPx3wotz504BAFxc1A+seh9TUykA4NWrl3l6PBU+LLZEpLGLF/+GEAKWllY4d+405HI5kpOf4sKFs9i2bT0SE+MxY8ZS9OrVP8d5WFpaoWbNOh+c5fDh/XjxIg1paakA3pyAYu/eN8fJNmvWBmZm5jh16hi6dWuGkSMnYeTISQAAmUwGf/8mKvPbtm09jIyMlab16zcCffp0QPfuLdCv30jY2dnj779PY8mScFSsWBUBAa3zlL1y5eoAgNWrv0OXLsEwMTFBuXKVYGlplaf5UcHHYktEGss+c9SaNYuwZs0iSKVS2NjYoly5yujaNQSfffYlnJxc8iXLhAkDce/eXcX9vXt/xt69PwP4b/dz9mFGcrnqWbA0ERj4CbZuPYTvv5+JyZOHIyUlGa6uZdGzZ38MGTIh1+Nsc+Pv3wRDhkzA9u0bsGnTKsjlcvz88xG1XwKoaJAInmlcSUpKCmxsbHD1ajKsrKwNHYcoX/1x4hK695ipuH/5wgqVQ3+IioPU1BRUrmyD5ORkWFt/eC3gGaSIiIj0jMWWiIhIz/ibLVExcP9+LJ48efzefjG3bwLiieL+lSvnYWVlptFz2NnZK86MRETKWGyJirj792PRuHEVja+/+rYuXfZr3NfMzBzHjv3LgkukBostURH35MljvHz5AqGhofDw8NDLc9y5cwehoaF48uQxiy2RGiy2RMWEh4cHKleubOgYRMUSB0gRERHpGYstERGRnhWqYnv8+HEEBQXB1dUVEokEu3btUpouhEBoaChcXV1hZmaGJk2a4PLly4YJS0RE9P8KVbFNS0tDzZo1sWTJErXTZ8+ejfnz52PJkiWIioqCs7MzWrRogdTU1HxOSkRE9J9CNUCqdevWaN1a/Ym/hRBYuHAhvvnmG3Tq1AkAsGHDBjg5OWHTpk3o3z/nE6MTERHpU6Hass1NTEwMEhISEBgYqGiTSqVo3LgxTp48mePj0tPTkZKSonQjIiLSpSJTbBMSEgAATk5OSu1OTk6KaeqEh4fDxsZGcStbtqxecxIRUfFTZIptNolEonRfCKHS9rYJEyYgOTlZcYuLi9N3RCIiKmYK1W+2uXF2dgbwZgvXxeW/62kmJiaqbO2+TSqVQiqV6j0fEREVX0Vmy9bT0xPOzs6IjIxUtGVkZODYsWPw9/c3YDIiIiruCtWW7fPnz3Hz5k3F/ZiYGJw/fx52dnZwc3PDiBEjMGPGDFSoUAEVKlTAjBkzYG5ujh49ehgwNRERFXeFqtiePXsWAQEBivujRo0CAAQHB2P9+vUYN24cXr58iUGDBuHp06eoV68eDhw4ACsrK0NFJiIiKlzFtkmTJhBC5DhdIpEgNDQUoaGh+ReKiIjoPYrMb7ZEREQFlUZbthEREVrPuEWLFjAzM9P6cUREREWNRsW2Q4cOWs1UIpHgxo0b8PLyyksmIiKiIkXj3cgJCQmQy+Ua3czNzfWZmYiIqFDRqNgGBwdrtUu4Z8+esLa2znMoIiKiokSj3cjr1q3TaqbLli3LUxgiIqKiSOvRyH369FF7fdi0tDT06dNHJ6GIiIiKEq2L7YYNG/Dy5UuV9pcvX+KHH37QSSgiIqKiROOTWqSkpEAIASEEUlNTIZPJFNOysrLw66+/wtHRUS8hiYiICjONi23JkiUhkUggkUhQsWJFlekSiQRhYWE6DUdERFQUaFxsjxw5AiEEmjZtil9++QV2dnaKaaampnB3d4erq6teQhIRERVmGhVbOzs7XL9+Hfb29ggODkbz5s15cn8iIiINaTRAKiMjAykpKQCAH374Aa9evdJrKCIioqJEoy1bPz8/dOjQAb6+vhBCYNiwYTme5GLt2rU6DUhERFTYaVRsN27ciAULFuDWrVuQSCRITk7m1i0REZGGNCq2Tk5OmDlzJgDA09MTP/74I0qVKqXXYEREREWF1hePj4mJ0UcOIiKiIkujAVKLFi3Sarfx8uXL1Z7SkYiIqDjSqNiOHDlSq+I5btw4PHr0KM+hiIiIihKNdiMLIdCsWTOUKKHZXmd1504mIiIqrjSqnpMnT9Zqpu3bt1c6wxQREVFxppdiS0RERP/R+hJ7REREpB0WWyIiIj1jsSUiItIzFlsiIiI907rYTpkyBS9evFBpf/nyJaZMmaKTUEREREWJ1qdrDAsLw4ABA2Bubq7U/uLFC4SFhWHSpEk6C2dIL16kwdjY2NAxiD7Yq1dvjntPT0/X2zHw6enpiud68SJNL89BlJ90/T6WCCGENg8wMjLCw4cP4eDgoNR++PBhdOvWrdCfOSolJQU2NjaGjkFERAVAcnIyrK2tP3g+Gm/Z2traQiKRQCKRoGLFipBIJIppWVlZeP78OQYMGPDBgYiIiIoajYvtwoULIYRAnz59EBYWprT1Z2pqCg8PD/j5+eklpCEc2fwtLM1lho5BpBPxj54iOeXNWIsXFjV1Nl/ztH8U/7exNoeLg63O5k1kSM9fvELAp1N1Nj+Ni21wcDCAN9ez9ff3h4mJic5CFETmMlOYm5kaOgaRTpRzc1L8P9Vad8XWKiVFZ/MiKkjkcrlO56f1AKnGjRtDLpfj+vXrSExMVAnUqFEjnYUjIiIqCrQutqdPn0aPHj1w9+5dvDu2SiKRICsrS2fhiIiIigKtj7MdMGAA6tSpg0uXLuHJkyd4+vSp4vbkyRN9ZNTa0qVL4enpCZlMBl9fX/zxxx+GjkRERMWY1lu2N27cwPbt21G+fHl95PlgW7duxYgRI7B06VI0aNAAK1asQOvWrXHlyhW4ubkZOh4RERVDWhfbevXq4ebNmwW22M6fPx99+/bFF198AeDNKOrff/8dy5YtQ3h4uMbzefHKBEZGRXsQGBVPL010d5ZW45f8G6Gi6cUr3f4kqlGxvXDhguL/Q4cOxejRo5GQkIDq1aurjEquUaOGTgNqIyMjA+fOncP48eOV2gMDA3Hy5Em1j0lPT1ec/QZ4c1ILAAj4dCCADz+Qmahoq2voAER6kgJg/Ht7aUqjYuvj4wOJRKI0IKpPnz6K/2dPM/QAqcePHyMrKwtOTk5K7U5OTkhISFD7mPDwcISFheVHPCIiKqY0KrYxMTH6zqFTb5/dCoDii4A6EyZMwKhRoxT3U1JSULZsWRzZvIwntaAi6W5mJY3P9jZv3jylsQ4HDx7EihUrFPeNRAZKlbTED/MH6zwnkSG9OamF7uanUbF1d3fX3TPqkb29PYyNjVW2YhMTE1W2drNJpVJIpVKVdnNZJszNeCECKnpMJZl4/PiuRn1NTDJgZvbfsfRyearKYyUSG5ibZeo0I5GhyeW6fU9rPUAqIiJCbbtEIoFMJkP58uXh6en5wcHywtTUFL6+voiMjETHjh0V7ZGRkWjfvr1BMhEVNEZGRioXEsnJu1e+kslkiscmJSXp/Cw7REWV1sW2Q4cOKr/fAsq/2zZs2BC7du2CrW3+nyd11KhR6NWrF+rUqQM/Pz+sXLkSsbGxvEgC0f+zt7fHnj178vTYli1bomXLlgCAoKCgQn+VL6L8ovUxAJGRkahbty4iIyORnJyM5ORkREZG4qOPPsLevXtx/PhxJCUlYcyYMfrI+17dunXDwoULMWXKFPj4+OD48eP49ddfC82ucCIiKnq03rIdPnw4Vq5cCX9/f0Vbs2bNIJPJ8OWXX+Ly5ctYuHCh0mjl/DZo0CAMGjTIYM9PVFBt23caT0UczM3N0aFDB0PHISo2tC62t27dUnshXWtra9y+fRsAUKFCBTx+/PjD0xGRTi3/6SAePk6Gg4MDiy1RPtJ6N7Kvry/Gjh2r9FvNo0ePMG7cONSt++YA9xs3bqBMmTK6S0lEBU6pUqXg6GAHe1srQ0chKvC03rJds2YN2rdvjzJlyqBs2bKQSCSIjY2Fl5cXdu/eDQB4/vw5vv32W52HJaKCY/369TBNvwtpeqyhoxAVeFoX20qVKuHff//F77//juvXr0MIgcqVK6NFixYwMnqzoczdU0RERP/RutgCbw7zadWqFVq1aqXrPEREREWORsV20aJF+PLLLyGTybBo0aJc+w4bNkwnwYiIiIoKjYrtggUL8Nlnn0Emk2HBggU59pNIJCy2RMXEzJkz8fxpPGwtgNARnQ0dh6hA0/pCBIXtogREpB9//vknHj16BCd7G0NHISrw8nwV6YyMDFy7dg2vX7/WZR4iIqIiR+ti++LFC/Tt2xfm5uaoVq0aYmPfDPsfNmwYZs6cqfOARKQ77qXt4enpqXTZPCLSP62L7YQJE/DPP//g6NGjkMn+u95r8+bNsXXrVp2GIyLdWjdnADZv3ozvv//e0FGIihWtD/3ZtWsXtm7divr16ytdkL1q1aq4deuWTsMREREVBVpv2T569AiOjo4q7WlpaUrFl4iIiN7QutjWrVsX+/btU9zPLrCrVq2Cn5+f7pIREREVEVrvRg4PD0erVq1w5coVvH79Gt999x0uX76MU6dO4dixY/rISEQ6Mi58Ex6/2A4bGxtMmTLF0HGIig2tt2z9/f3x559/4sWLFyhXrhwOHDgAJycnnDp1Cr6+vvrISEQ6cvbibZw5cwbR0dGGjkJUrOTp3MjVq1fHhg0bdJ2FiAqRwMBApD17AFszuaGjEBV4eSq2crkcN2/eRGJiIuRy5T+0Ro0a6SQYERVsQ4cO5SX2iDSkdbE9ffo0evTogbt370IIoTRNIpEgKytLZ+GIiIiKAq2L7YABA1CnTh3s27cPLi4uPNyHiIjoPbQutjdu3MD27dtRvnx5feQhIiIqcrQutvXq1cPNmzdZbImKuW7duuHxo0Q42lli79pxho5DVKBpVGwvXLig+P/QoUMxevRoJCQkoHr16jAxMVHqW6NGDd0mJKIC6cWLF0h78RIvzE0NHYWowNOo2Pr4+EAikSgNiOrTp4/i/9nTOECKiIhIldYXjyeiwqtz64+QlFkKlpaWho5CVKxoVGzd3d31nYOI8sGgXoFItf7Y0DGIih2tT9dIRERE2mGxJSIi0jMWWyIiIj3L07mRiahwatpjGh4+ToaDgwP27Nlj6DhExUaetmyfPXuG1atXY8KECXjy5AkA4O+//8b9+/d1Go6IiKgo0HrL9sKFC2jevDlsbGxw584d9OvXD3Z2dti5cyfu3r2LH374QR85iaiA+eqrr5D1/B4sjZINHYWowNN6y3bUqFEICQnBjRs3IJPJFO2tW7fG8ePHdRqOiAquhg0bokWAH5rUr2roKEQFntbFNioqCv3791dpL126NBISEnQSSp3p06fD398f5ubmKFmypNo+sbGxCAoKgoWFBezt7TFs2DBkZGToLRMREZEmtN6NLJPJkJKSotJ+7do1ODg46CSUOhkZGejSpQv8/PywZs0alelZWVlo27YtHBwccOLECSQlJSE4OBhCCCxevFhvuYiIiN5H62Lbvn17TJkyBdu2bQPw5rzIsbGxGD9+PP73v//pPGC2sLAwAMD69evVTj9w4ACuXLmCuLg4uLq6AgDmzZuHkJAQTJ8+HdbW1nrLRlQcXb16FSItDhbiMapVLGPoOEQFmta7kefOnYtHjx7B0dERL1++ROPGjVG+fHlYWVlh+vTp+siokVOnTsHb21tRaAGgZcuWSE9Px7lz53J8XHp6OlJSUpRuRPR+Y8eORe/B32Jo6HpDRyEq8LTesrW2tsaJEydw+PBh/P3335DL5ahduzaaN2+uj3waS0hIgJOTk1Kbra0tTE1Nc/0tOTw8XLHVTEREpA9ab9neuXMHANC0aVOMGTMG48aNy3OhDQ0NhUQiyfV29uxZjecnkUhU2rIv/ZeTCRMmIDk5WXGLi4vL07IQERHlROstWy8vL/j7+6NXr17o0qUL7Ozs8vzkQ4YMQffu3XPt4+HhodG8nJ2dcebMGaW2p0+fIjMzU2WL921SqRRSqVSj5yAq7GZ+9SmSTSrDxMTE0FGIihWti+3Zs2exefNmTJs2DcOHD0fLli3Rs2dPfPLJJ1oXLXt7e9jb22sbQS0/Pz9Mnz4d8fHxcHFxAfBm0JRUKoWvr69OnoOosPuoZjmkWtc3dAyiYkfr3ci1a9fGnDlzEBsbi/3798PR0RH9+/eHo6Mj+vTpo4+MAN4cQ3v+/HnExsYiKysL58+fx/nz5/H8+XMAQGBgIKpWrYpevXohOjoahw4dwpgxY9CvXz+ORCYiIoPK81V/JBIJAgICsGrVKhw8eBBeXl7YsGGDLrMpmTRpEmrVqoXJkyfj+fPnqFWrFmrVqqX4TdfY2Bj79u2DTCZDgwYN0LVrV3To0AFz587VWyYiIiJN5PmqP3Fxcdi8eTM2bdqEixcvws/PD0uWLNFlNiXr16/P8RjbbG5ubti7d6/eMhAVdn/9cwvJJiYwMTHhzytE+UjrYrty5Ur89NNP+PPPP1GpUiV89tln2LVrl8YDmYjIcMbP2sxL7BEZgNbFdurUqejevTu+++47+Pj46CESERFR0aJ1sY2Njc31uFUiKh62bNkCk/RYyDLuGToKUYGnUbG9cOECvL29YWRkhIsXL+bat0aNGjoJRkQFm4WFBUxLmENaQvb+zkTFnEbF1sfHBwkJCXB0dISPjw8kEgmEEIrp2fclEgmysrL0FpaIiKgw0qjYxsTEKC6fFxMTo9dARERERY1Gxdbd3V3x/7t378Lf3x8lSig/9PXr1zh58qRSXyIqujZt2oRXKQ9gY/oKIZ0bGzoOUYGm9QCpgIAAxMfHw9HRUak9OTkZAQEB3I1MVExs3rwZjx49gpO9DYst0XtofQapnK6ik5SUBAsLC52EIiIiKko03rLt1KkTgDeDoUJCQpQuOpCVlYULFy7A399f9wmJiIgKOY2LrY2NDYA3W7ZWVlYwMzNTTDM1NUX9+vXRr18/3SckIp05vGkiUq0/NnQMomJH42K7bt06AG+uLztmzBjuMiYqpDZt2oTNmze/t1+lSpVULuQxZswYXLt2DcCbn46ISDNaD5CaPHmyPnIQUT7JeHYbjx49em8/Z4eSME2/q9SW/OShRo8lImV5uurP9u3bsW3bNsTGxiIjI0Np2t9//62TYESkH1YywMne5r39SlmbQJoeq9L27mPtba10mo+oKNK62C5atAjffPMNgoODsXv3bvTu3Ru3bt1CVFQUBg8erI+MRKRDIZ0b5/lQne+n9NZxGqLiQetDf5YuXYqVK1diyZIlMDU1xbhx4xAZGYlhw4YhOTlZHxmJiIgKNa2LbWxsrOIQHzMzM6SmpgIAevXqpdGgCyIiouJG62Lr7OysGIXo7u6O06dPA3hzzuS3L05AREREb2hdbJs2bYo9e/YAAPr27YuRI0eiRYsW6NatGzp27KjzgERERIWd1gOkVq5cCblcDgAYMGAA7OzscOLECQQFBWHAgAE6D0hERFTYaV1sjYyMYGT03wZx165d0bVrV52GIiIiKko0KrYXLlzQeIY1atTIcxgiIqKiSKNi6+PjA4lE8t4BUBKJhJfYIyIieodGxTYmJkbfOYiIiIosjYqtu7u7vnMQEREVWVof+gMAP/74Ixo0aABXV1fcvfvmROULFy7E7t27dRqOiIioKNC62C5btgyjRo1CmzZt8OzZM8VvtCVLlsTChQt1nY+IiKjQ07rYLl68GKtWrcI333wDY2NjRXudOnVw8eJFnYYjIiIqCrQutjExMahVq5ZKu1QqRVpamk5CERERFSVaF1tPT0+cP39epX3//v2oWrWqLjIREREVKVqfQWrs2LEYPHgwXr16BSEE/vrrL2zevBnh4eFYvXq1PjISEREValoX2969e+P169cYN24cXrx4gR49eqB06dL47rvv0L17d31kJCIiKtS0LrYA0K9fP/Tr1w+PHz+GXC6Ho6MjAOD+/fsoXbq0TgMSEREVdnk6zjabvb09HB0dkZCQgKFDh6J8+fK6ykVERFRkaFxsnz17hs8++wwODg5wdXXFokWLIJfLMWnSJHh5eeH06dNYu3atXkLeuXMHffv2haenJ8zMzFCuXDlMnjwZGRkZSv1iY2MRFBQECwsL2NvbY9iwYSp9iIiI8pvGu5G//vprHD9+HMHBwfjtt98wcuRI/Pbbb3j16hX279+Pxo0b6y3k1atXIZfLsWLFCpQvXx6XLl1Cv379kJaWhrlz5wIAsrKy0LZtWzg4OODEiRNISkpCcHAwhBBYvHix3rIRERG9j0S871I+/8/d3R1r1qxB8+bNcfv2bZQvXx7Dhg0z2Fmj5syZg2XLluH27dsA3hx61K5dO8TFxcHV1RUAsGXLFoSEhCAxMRHW1tZq55Oeno709HTF/ZSUFJQtWxZndk6FpYVM/wtCREQFzvO0V6jX8VskJyfnWD+0ofFu5AcPHiiOo/Xy8oJMJsMXX3zxwQHyKjk5GXZ2dor7p06dgre3t6LQAkDLli2Rnp6Oc+fO5Tif8PBw2NjYKG5ly5bVa24iIip+NC62crkcJiYmivvGxsawsLDQS6j3uXXrFhYvXowBAwYo2hISEuDk5KTUz9bWFqampkhISMhxXhMmTEBycrLiFhcXp7fcRERUPGn8m60QAiEhIZBKpQCAV69eYcCAASoFd8eOHRo/eWhoKMLCwnLtExUVhTp16ijuP3jwAK1atUKXLl1UtqwlEona3Oras0mlUsUyERER6YPGxTY4OFjpfs+ePT/4yYcMGfLeE2F4eHgo/v/gwQMEBATAz88PK1euVOrn7OyMM2fOKLU9ffoUmZmZKlu8RERE+UnjYrtu3TqdP7m9vT3s7e016nv//n0EBATA19cX69atg5GR8h5wPz8/TJ8+HfHx8XBxcQEAHDhwAFKpFL6+vjrPTkREpKk8nUEqvz148ABNmjSBm5sb5s6di0ePHimmOTs7AwACAwNRtWpV9OrVC3PmzMGTJ08wZswY9OvXTycjyYiIiPKqUBTbAwcO4ObNm7h58ybKlCmjNC37yCVjY2Ps27cPgwYNQoMGDWBmZoYePXoojsMlIiIyFI2Psy0uUlJSYGNjw+NsiYiKMYMdZ0tERER5w2JLRESkZyy2REREesZiS0REpGcstkRERHrGYktERKRnLLZERER6xmJLRESkZyy2REREesZiS0REpGcstkRERHrGYktERKRnLLZERER6xmJLRESkZyy2REREesZiS0REpGcstkRERHrGYktERKRnLLZERER6xmJLRESkZyy2REREesZiS0REpGcstkRERHrGYktERKRnLLZERER6xmJLRESkZyy2REREesZiS0REpGcstkRERHrGYktERKRnLLZERER6xmJLRESkZyy2REREelZoiu0nn3wCNzc3yGQyuLi4oFevXnjw4IFSn9jYWAQFBcHCwgL29vYYNmwYMjIyDJSYiIjojUJTbAMCArBt2zZcu3YNv/zyC27duoXOnTsrpmdlZaFt27ZIS0vDiRMnsGXLFvzyyy8YPXq0AVMTEREBEiGEMHSIvIiIiECHDh2Qnp4OExMT7N+/H+3atUNcXBxcXV0BAFu2bEFISAgSExNhbW2t0XxTUlJgY2ODMzunwtJCps9FICKiAup52ivU6/gtkpOTNa4fuSk0W7Zve/LkCX766Sf4+/vDxMQEAHDq1Cl4e3srCi0AtGzZEunp6Th37lyO80pPT0dKSorSjYiISJcKVbH96quvYGFhgVKlSiE2Nha7d+9WTEtISICTk5NSf1tbW5iamiIhISHHeYaHh8PGxkZxK1u2rN7yExFR8WTQYhsaGgqJRJLr7ezZs4r+Y8eORXR0NA4cOABjY2N8/vnneHsvuEQiUXkOIYTa9mwTJkxAcnKy4hYXF6fbhSQiomKvhCGffMiQIejevXuufTw8PBT/t7e3h729PSpWrIgqVaqgbNmyOH36NPz8/ODs7IwzZ84oPfbp06fIzMxU2eJ9m1QqhVQq/aDlICIiyo1Bi2128cyL7C3a9PR0AICfnx+mT5+O+Ph4uLi4AAAOHDgAqVQKX19f3QQmIiLKA4MWW0399ddf+Ouvv9CwYUPY2tri9u3bmDRpEsqVKwc/Pz8AQGBgIKpWrYpevXphzpw5ePLkCcaMGYN+/frpZCQZERFRXhWKAVJmZmbYsWMHmjVrhkqVKqFPnz7w9vbGsWPHFLuAjY2NsW/fPshkMjRo0ABdu3ZFhw4dMHfuXAOnJyKi4q5QbNlWr14dhw8ffm8/Nzc37N27Nx8SERERaa5QbNkSEREVZoViyzY/ZQ+8ev7ilYGTEBGRoWTXAF2dZLHQnq5RX27fvo1y5coZOgYRERUAt27dgpeX1wfPh1u277CzswPw5gpCNjY2Bk6jvZSUFJQtWxZxcXGFdhR2YV8G5je8wr4MhT0/UPiXITk5GW5uboqa8KFYbN9hZPTmZ2wbG5tC+QbJZm1tXajzA4V/GZjf8Ar7MhT2/EDhX4bsmvDB89HJXIiIiChHLLZERER6xmL7DqlUismTJxfa8yUX9vxA4V8G5je8wr4MhT0/UPiXQdf5ORqZiIhIz7hlS0REpGcstkRERHrGYktERKRnLLZERER6xmIL4M6dO+jbty88PT1hZmaGcuXKYfLkycjIyFDqFxsbi6CgIFhYWMDe3h7Dhg1T6WNI06dPh7+/P8zNzVGyZEm1fSQSicpt+fLl+Rs0B5rkL+jr4F0eHh4qr/f48eMNHStXS5cuhaenJ2QyGXx9ffHHH38YOpJGQkNDVV5rZ2dnQ8fK1fHjxxEUFARXV1dIJBLs2rVLaboQAqGhoXB1dYWZmRmaNGmCy5cvGyasGu/LHxISorJO6tevb5iwaoSHh6Nu3bqwsrKCo6MjOnTogGvXrin10dU6YLEFcPXqVcjlcqxYsQKXL1/GggULsHz5cnz99deKPllZWWjbti3S0tJw4sQJbNmyBb/88gtGjx5twOTKMjIy0KVLFwwcODDXfuvWrUN8fLziFhwcnE8Jc/e+/IVhHagzZcoUpdd74sSJho6Uo61bt2LEiBH45ptvEB0djY8//hitW7dGbGysoaNppFq1akqv9cWLFw0dKVdpaWmoWbMmlixZonb67NmzMX/+fCxZsgRRUVFwdnZGixYtkJqams9J1XtffgBo1aqV0jr59ddf8zFh7o4dO4bBgwfj9OnTiIyMxOvXrxEYGIi0tDRFH52tA0FqzZ49W3h6eiru//rrr8LIyEjcv39f0bZ582YhlUpFcnKyISLmaN26dcLGxkbtNABi586d+ZpHWznlL0zrIJu7u7tYsGCBoWNo7KOPPhIDBgxQaqtcubIYP368gRJpbvLkyaJmzZqGjpFn7/5tyuVy4ezsLGbOnKloe/XqlbCxsRHLly83QMLcqftsCQ4OFu3btzdInrxITEwUAMSxY8eEELpdB9yyzUFycrLSCahPnToFb29vuLq6KtpatmyJ9PR0nDt3zhAR82zIkCGwt7dH3bp1sXz5csjlckNH0khhXQezZs1CqVKl4OPjg+nTpxfY3d4ZGRk4d+4cAgMDldoDAwNx8uRJA6XSzo0bN+Dq6gpPT090794dt2/fNnSkPIuJiUFCQoLS+pBKpWjcuHGhWR8AcPToUTg6OqJixYro168fEhMTDR0pR8nJyQD+uyCNLtcBL0Sgxq1bt7B48WLMmzdP0ZaQkAAnJyelfra2tjA1NUVCQkJ+R8yzqVOnolmzZjAzM8OhQ4cwevRoPH78uEDv2sxWGNfB8OHDUbt2bdja2uKvv/7ChAkTEBMTg9WrVxs6morHjx8jKytL5TV2cnIqsK/v2+rVq4cffvgBFStWxMOHDzFt2jT4+/vj8uXLKFWqlKHjaS37NVe3Pu7evWuISFpr3bo1unTpAnd3d8TExODbb79F06ZNce7cuQJ3ZikhBEaNGoWGDRvC29sbgG7XQZHeslU3YOLd29mzZ5Ue8+DBA7Rq1QpdunTBF198oTRNIpGoPIcQQm27IZchNxMnToSfnx98fHwwevRoTJkyBXPmzCk0+Q2xDt6lzTKNHDkSjRs3Ro0aNfDFF19g+fLlWLNmDZKSkvItr7befS3z+/XNq9atW+N///sfqlevjubNm2Pfvn0AgA0bNhg42YcprOsDALp164a2bdvC29sbQUFB2L9/P65fv65YNwXJkCFDcOHCBWzevFllmi7WQZHesh0yZAi6d++eax8PDw/F/x88eICAgAD4+flh5cqVSv2cnZ1x5swZpbanT58iMzNT5VuPLmm7DNqqX78+UlJS8PDhQ70shy7zG2odvOtDlil7JObNmzcL3NaWvb09jI2NVbZiExMT8/X11RULCwtUr14dN27cMHSUPMkeSZ2QkAAXFxdFe2FdHwDg4uICd3f3ArdOhg4dioiICBw/fhxlypRRtOtyHRTpYmtvbw97e3uN+t6/fx8BAQHw9fXFunXrVK5h6Ofnh+nTpyM+Pl7xoh84cABSqRS+vr46z55Nm2XIi+joaMhkshwPtflQusxvqHXwrg9ZpujoaABQ+sMtKExNTeHr64vIyEh07NhR0R4ZGYn27dsbMFnepKen499//8XHH39s6Ch54unpCWdnZ0RGRqJWrVoA3vyufuzYMcyaNcvA6fImKSkJcXFxBeb9L4TA0KFDsXPnThw9ehSenp5K03W6DnQzhqtwu3//vihfvrxo2rSpuHfvnoiPj1fcsr1+/Vp4e3uLZs2aib///lscPHhQlClTRgwZMsSAyZXdvXtXREdHi7CwMGFpaSmio6NFdHS0SE1NFUIIERERIVauXCkuXrwobt68KVatWiWsra3FsGHDDJz8jfflLwzr4G0nT54U8+fPF9HR0eL27dti69atwtXVVXzyySeGjpajLVu2CBMTE7FmzRpx5coVMWLECGFhYSHu3Llj6GjvNXr0aHH06FFx+/Ztcfr0adGuXTthZWVVoLOnpqYq3ucAFO+Xu3fvCiGEmDlzprCxsRE7duwQFy9eFJ9++qlwcXERKSkpBk7+Rm75U1NTxejRo8XJkydFTEyMOHLkiPDz8xOlS5cuMPkHDhwobGxsxNGjR5U+91+8eKHoo6t1wGIr3hxqAkDt7W13794Vbdu2FWZmZsLOzk4MGTJEvHr1ykCpVQUHB6tdhiNHjgghhNi/f7/w8fERlpaWwtzcXHh7e4uFCxeKzMxMwwb/f+/LL0TBXwdvO3funKhXr56wsbERMplMVKpUSUyePFmkpaUZOlquvv/+e+Hu7i5MTU1F7dq1FYdBFHTdunUTLi4uwsTERLi6uopOnTqJy5cvGzpWro4cOaL2PR8cHCyEeHPoyeTJk4Wzs7OQSqWiUaNG4uLFi4YN/Zbc8r948UIEBgYKBwcHYWJiItzc3ERwcLCIjY01dGyFnD73161bp+ijq3XAS+wRERHpWZEejUxERFQQsNgSERHpGYstERGRnrHYEhER6RmLLRERkZ6x2BIREekZiy0REZGesdgSERHpGYstUQElkUiwa9cuQ8fQuZCQEMXVkd5evqtXr6J+/fqQyWTw8fHJ8fFNmjRRPP78+fN6z0ukCyy2RPno7UJjYmICJycntGjRAmvXroVcLlfqGx8fj9atW2s038JWmFu1aqWyfJMnT4aFhQWuXbuGQ4cOYf369WovkLFjxw789ddf+ZiW6MOx2BLls+xCc+fOHezfvx8BAQEYPnw42rVrh9evXyv6OTs7F7gLbOuKVCpVWb5bt26hYcOGcHd3z/Xyg3Z2dnBwcMiPmEQ6w2JLlM+yC03p0qVRu3ZtfP3119i9ezf279+P9evXK/q9vbWakZGBIUOGwMXFBTKZDB4eHggPDwfw37VzO3bsCIlEorh/69YttG/fHk5OTrC0tETdunVx8OBBpSweHh6YMWMG+vTpAysrK7i5ualcy/nevXvo3r077OzsYGFhgTp16ihdV3jPnj3w9fWFTCaDl5cXwsLClL40aEIikeDcuXOYMmUKJBIJmjRpgt69eyM5OVmxJyA0NFSreRIVJCy2RAVA06ZNUbNmTezYsUPt9EWLFiEiIgLbtm3DtWvXsHHjRkVRjYqKAgCsW7cO8fHxivvPnz9HmzZtcPDgQURHR6Nly5YICgpCbGys0rznzZuHOnXqIDo6GoMGDcLAgQNx9epVxTwaN26MBw8eICIiAv/88w/GjRun2OX9+++/o2fPnhg2bBiuXLmCFStWYP369Zg+fbpWyx8fH49q1aph9OjRiI+PR0REBBYuXAhra2vEx8cjPj4eY8aM0WqeRAVJkb54PFFhUrlyZVy4cEHttNjYWFSoUAENGzaERCKBu7u7Ylr2LtWSJUvC2dlZ0V6zZk3UrFlTcX/atGnYuXMnIiIiMGTIEEV7mzZtMGjQIADAV199hQULFuDo0aOoXLkyNm3ahEePHiEqKgp2dnYAgPLlyyseO336dIwfPx7BwcEAAC8vL0ydOhXjxo3D5MmTNV52Z2dnlChRApaWloplsLGxgUQiUVomosKKxZaogBBCQCKRqJ0WEhKCFi1aoFKlSmjVqhXatWuHwMDAXOeXlpaGsLAw7N27Fw8ePMDr16/x8uVLlS3bGjVqKP6fXdwSExMBAOfPn0etWrUUhfZd586dQ1RUlNKWbFZWFl69eoUXL17A3Nxco2UnKupYbIkKiH///Reenp5qp9WuXRsxMTHYv38/Dh48iK5du6J58+bYvn17jvMbO3Ysfv/9d8ydOxfly5eHmZkZOnfujIyMDKV+JiYmSvclEoliN7GZmVmumeVyOcLCwtCpUyeVaTKZLNfHEhUnLLZEBcDhw4dx8eJFjBw5Msc+1tbW6NatG7p164bOnTujVatWePLkCezs7GBiYoKsrCyl/n/88QdCQkLQsWNHAG9+f71z545WuWrUqIHVq1crnuddtWvXxrVr15R2LeuKqampyjIRFVYstkT5LD09HQkJCcjKysLDhw/x22+/ITw8HO3atcPnn3+u9jELFiyAi4sLfHx8YGRkhJ9//hnOzs6K41A9PDxw6NAhNGjQAFKpFLa2tihfvjx27NiBoKAgSCQSfPvttyrH8r7Pp59+ihkzZqBDhw4IDw+Hi4sLoqOj4erqCj8/P0yaNAnt2rVD2bJl0aVLFxgZGeHChQu4ePEipk2b9kGvk4eHB54/f45Dhw6hZs2aMDc3525pKrQ4Gpkon/32229wcXGBh4cHWrVqhSNHjmDRokXYvXs3jI2N1T7G0tISs2bNQp06dVC3bl3cuXMHv/76K4yM3vwJz5s3D5GRkShbtixq1aoF4E2BtrW1hb+/P4KCgtCyZUvUrl1bq6ympqY4cOAAHB0d0aZNG1SvXh0zZ85U5GzZsiX27t2LyMhI1K1bF/Xr18f8+fOVBnDllb+/PwYMGIBu3brBwcEBs2fP/uB5EhmKRAghDB2CiIqPkJAQPHv27IPOeHXnzh14enoiOjo611M7EhUU3LIlony3d+9eWFpaYu/evVo/tnXr1qhWrZoeUhHpD7dsiShfJSYmIiUlBQDg4uICCwsLrR5///59vHz5EgDg5uYGU1NTnWck0jUWWyIiIj3jbmQiIiI9Y7ElIiLSMxZbIiIiPWOxJSIi0jMWWyIiIj1jsSUiItIzFlsiIiI9Y7ElIiLSs/8Dw45VgOUuU6wAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "\n", "##Now printing the conceptual model figure:\n", "\n", "fig = plt.figure()\n", "ax = fig.add_subplot(1, 1, 1)\n", "# sky\n", "sky = plt.Rectangle((-20, 2), width=50, height=20, fc=\"b\", zorder=0, alpha=0.1)\n", "ax.add_patch(sky)\n", "\n", "# Aquifer:\n", "ground = plt.Rectangle(\n", " (-20, -32.57),\n", " width=50,\n", " height=34.57,\n", " fc=np.array([209, 179, 127]) / 255,\n", " zorder=0,\n", " alpha=0.9,\n", ")\n", "ax.add_patch(ground)\n", "\n", "well = plt.Rectangle(\n", " (-1, -(0.47 + 13.8)),\n", " width=2,\n", " height=(0.47 + 13.8) + 2,\n", " fc=np.array([200, 200, 200]) / 255,\n", " zorder=1,\n", ")\n", "ax.add_patch(well)\n", "\n", "# Wellhead\n", "wellhead = plt.Rectangle(\n", " (-1.25, 2),\n", " width=2.5,\n", " height=10,\n", " fc=np.array([200, 200, 200]) / 255,\n", " zorder=2,\n", " ec=\"k\",\n", ")\n", "ax.add_patch(wellhead)\n", "\n", "# Screen for the well:\n", "screen = plt.Rectangle(\n", " (-1, -(0.47 + 13.8)),\n", " width=2,\n", " height=13.8,\n", " fc=np.array([200, 200, 200]) / 255,\n", " alpha=1,\n", " zorder=2,\n", " ec=\"k\",\n", " ls=\"--\",\n", ")\n", "screen.set_linewidth(2)\n", "ax.add_patch(screen)\n", "pumping_arrow = plt.Arrow(x=0, y=10, dx=0, dy=6, color=\"#00035b\")\n", "ax.add_patch(pumping_arrow)\n", "ax.text(x=0.5, y=13, s=r\"$ D = 1.48$ ft\", fontsize=\"large\")\n", "\n", "# last line\n", "line = plt.Line2D(xdata=[-200, 1200], ydata=[2, 2], color=\"k\")\n", "ax.add_line(line)\n", "\n", "# water table\n", "line = plt.Line2D(xdata=[-200, 1200], ydata=[0, 0], color=\"blue\")\n", "ax.add_line(line)\n", "\n", "\n", "ax.set_xlim([-20, 20])\n", "ax.set_ylim([-32, 20])\n", "ax.set_xlabel(\"Distance [ft]\")\n", "ax.set_ylabel(\"Relative height [ft]\")\n", "ax.set_title(\"Conceptual Model - Falling Head Example\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 2. Set basic parameters\n", "\n", "Parameters here declared are already converted from feet and inches to meters" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "rw = 0.127 # well radius, m\n", "rc = 0.0508 # well casing radius, m\n", "L = 4.20624 # screen length, m\n", "b = -9.9274 # aquifer thickness, m\n", "zt = -0.1433 # depth to top of the screen, m\n", "H0 = 0.4511 # initial displacement in the well, m\n", "zb = zt - L # bottom of the screen, m" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 3. Converting slug displacement to volume" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Slug: 0.00366 m^3\n" ] } ], "source": [ "Q = np.pi * rc**2 * H0\n", "print(f\"Slug: {Q:.5f} m^3\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 4. Load data\n", "\n", "Drawdown data is available in feet and seconds and are converted to meters and days" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "data = np.loadtxt(\"data/falling_head.txt\", skiprows=2)\n", "to = data[:, 0] / 60 / 60 / 24 # convert time from seconds to days\n", "ho = (10 - data[:, 1]) * 0.3048 # convert drawdown from ft to meters" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 5. Create First Model - three layers\n", "\n", "We begin with a model with just three layers. We arranged the layers to match the screen length. The first layer is located just above the screen, the second layer is located at the screen depths, and the last layer is just below the screen, up to the total aquifer depth.\n", "\n", "We set the model in the same manner as in [Slug 1 - Pratt County](slug1_pratt_county.ipynb)." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "self.neq 1\n", "solution complete\n" ] } ], "source": [ "ml_0 = ttm.Model3D(kaq=10, z=[0, zt, zb, b], Saq=1e-4, tmin=1e-5, tmax=0.01)\n", "w_0 = ttm.Well(ml_0, xw=0, yw=0, rw=rw, rc=rc, tsandQ=[(0, -Q)], layers=1, wbstype=\"slug\")\n", "ml_0.solve()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 6. Model calibration\n", "\n", "The procedures for calibration can be seen in [Unconfined 1 - Vennebulten](unconfined1_vennebulten.ipynb)\n", "\n", "We calibrate hydraulic conductivity and specific storage, as in the KGS model (Hyder et al. 1994)." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "...............................................................\n", "Fit succeeded.\n", "[[Fit Statistics]]\n", " # fitting method = leastsq\n", " # function evals = 60\n", " # data points = 27\n", " # variables = 2\n", " chi-square = 0.00114578\n", " reduced chi-square = 4.5831e-05\n", " Akaike info crit = -267.822650\n", " Bayesian info crit = -265.230976\n", "[[Variables]]\n", " kaq0_2: 0.59981087 +/- 0.03385855 (5.64%) (init = 10)\n", " Saq0_2: 2.0901e-04 +/- 6.5392e-05 (31.29%) (init = 0.0001)\n", "[[Correlations]] (unreported correlations are < 0.100)\n", " C(kaq0_2, Saq0_2) = -0.9720\n" ] } ], "source": [ "ca_0 = ttm.Calibrate(ml_0)\n", "ca_0.set_parameter(name=\"kaq0_2\", initial=10)\n", "ca_0.set_parameter(name=\"Saq0_2\", initial=1e-4, pmin=1e-7)\n", "ca_0.seriesinwell(name=\"obs\", element=w_0, t=to, h=ho)\n", "ca_0.fit(report=True)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
optimalstdperc_stdpminpmaxinitialparray
kaq0_20.5998110.0338595.644871-infinf10.0000[0.5998108693035121, 0.5998108693035121, 0.599...
Saq0_20.0002090.00006531.2868011.000000e-07inf0.0001[0.00020900772948173607, 0.0002090077294817360...
\n", "
" ], "text/plain": [ " optimal std perc_std pmin pmax initial \\\n", "kaq0_2 0.599811 0.033859 5.644871 -inf inf 10.0000 \n", "Saq0_2 0.000209 0.000065 31.286801 1.000000e-07 inf 0.0001 \n", "\n", " parray \n", "kaq0_2 [0.5998108693035121, 0.5998108693035121, 0.599... \n", "Saq0_2 [0.00020900772948173607, 0.0002090077294817360... " ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "RMSE: 0.006514317872522675\n" ] } ], "source": [ "display(ca_0.parameters)\n", "print(\"RMSE:\", ca_0.rmse())" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcoAAAE/CAYAAADGw4N2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABYiklEQVR4nO3deXhM59vA8e9ksosksoqIJKit9n2pWlpbqKJIKbG3tpYqSv1qq4pqa3m1aEtrqSpVVBUV+97aWqV2IUqIWBIEmWSe949pRkciZmKSyXJ/rutczJnnnHOfmZPcec55Fo1SSiGEEEKIDNnZOgAhhBAiN5NEKYQQQmRCEqUQQgiRCUmUQgghRCYkUQohhBCZkEQphBBCZEISpRBCCJEJSZRCCCFEJiRRCiGEEJmQRJlPLFiwAI1Gg0ajYdu2beneV0pRunRpNBoNjRs3tuqxNRoN48ePt3i78+fPo9FoWLBggVXjyW4Zxb1nzx7Gjx/PrVu3bBJTZscPCQmhTZs2OR9UFmT1WhLmady4cZZ//kNCQujZs6dV48krJFHmM4ULF2b+/Pnp1m/fvp2zZ89SuHBhG0SV/+3Zs4cJEybYNFHa8vhC5GeSKPOZ8PBwfvzxRxITE03Wz58/n3r16lGiRAkbRWZdOp2OlJQUW4eR7927dw8ZDtpAKcW9e/dsHYawAUmU+UyXLl0AWLp0qXFdQkICP/74I717985wmxs3bjBw4EACAwNxdHSkZMmSjBkzhgcPHpiUS0xMpF+/fnh7e+Pm5kbLli05depUhvs8ffo0Xbt2xc/PDycnJ8qXL8/nn3+epXPatm0bGo2GxYsX88477xAYGIiTkxNnzpwBYNOmTbzwwgu4u7vj6upKgwYN2Lx5s8k+rl27xuuvv05QUBBOTk74+vrSoEEDNm3aZCzzuFtLT7pdNX78eEaMGAFAaGhoulvgW7ZsoXHjxnh7e+Pi4kKJEiV45ZVXSEpKytLnYenx02zYsIHq1avj4uJCuXLl+Prrr03eT7t9v3HjRnr37o2vry+urq7G62DZsmXUq1ePQoUK4ebmRosWLTh8+HC6eA4cOEDbtm3x8vLC2dmZatWqsXz58iyd27Vr1xg4cCAVKlTAzc0NPz8/mjZtys6dO41llFI888wztGjRIt32d+7cwcPDg0GDBhnXJSYmMnz4cEJDQ3F0dCQwMJChQ4dy9+5dk201Gg2DBw9m7ty5lC9fHicnJxYuXAjAnDlzqFKlCm5ubhQuXJhy5crx3nvvZXouabfsP/74Yz766CNCQkJwcXGhcePGnDp1Cp1Ox6hRoyhWrBgeHh60b9+euLg4k33o9XqmTp1KuXLlcHJyws/Pj4iICP755x+Tckoppk6dSnBwMM7OzlSvXp3169dnGJe5n0eBpkS+8M033yhA7d+/X3Xv3l3Vrl3b+N6cOXNUoUKFVGJionr22WdVo0aNjO/du3dPVa5cWRUqVEh98sknauPGjer9999X9vb2KiwszFhOr9erJk2aKCcnJ/Xhhx+qjRs3qnHjxqmSJUsqQI0bN85Y9tixY8rDw0NVqlRJLVq0SG3cuFG98847ys7OTo0fP95YLjo6WgHqm2++yfTctm7dqgAVGBioOnbsqNasWaPWrl2rrl+/rhYvXqw0Go1q166dWrlypfr5559VmzZtlFarVZs2bTLuo0WLFsrX11d9+eWXatu2bWr16tVq7Nix6vvvvzeWCQ4OVj169Eh3/EaNGpl8Zo/GffHiRfXmm28qQK1cuVLt3btX7d27VyUkJKjo6Gjl7OysmjVrplavXq22bdumlixZorp3765u3ryZ6XmbK7Pjp51X8eLFVYUKFdSiRYvUr7/+qjp16qQAtX37duN+0q6hwMBA9frrr6v169erFStWqJSUFPXhhx8qjUajevfurdauXatWrlyp6tWrpwoVKqSOHTtm3MeWLVuUo6OjatiwoVq2bJnasGGD6tmzp1nfs1Iq3bV04sQJNWDAAPX999+rbdu2qbVr16o+ffooOzs7tXXrVmO5mTNnKo1Go06dOmWyv88//1wBxhjv3r2rqlatqnx8fNS0adPUpk2b1MyZM5WHh4dq2rSp0uv1JrEEBgaqypUrq++++05t2bJFHT16VC1dulQB6s0331QbN25UmzZtUnPnzlVvvfVWpueWdt0EBwerl156Sa1du1Z9++23yt/fX5UpU0Z1795d9e7dW61fv17NnTtXubm5qZdeeslkH6+//roC1ODBg9WGDRvU3Llzla+vrwoKClLXrl0zlhs3bpwCVJ8+fdT69evVl19+qQIDA1XRokVNrmVLPo/H/XwUBJIo84n/Jsq0xHL06FGllFK1atVSPXv2VEqpdIly7ty5ClDLly832d9HH32kALVx40allFLr169XgJo5c6ZJuQ8//DDdL7cWLVqo4sWLG39Rpxk8eLBydnZWN27cUEpZniiff/55k/V3795VXl5e6X6ZpKamqipVqpj8seDm5qaGDh2a6XGymiiVUurjjz9WgIqOjjbZdsWKFQpQf/zxR6bHflqPO75ShvNydnZWFy5cMK67d++e8vLyUm+88YZxXdo1FBERYbJ9TEyMsre3V2+++abJ+tu3b6uiRYuqzp07G9eVK1dOVatWTel0OpOybdq0UQEBASo1NTXT83j0WnpUSkqK0ul06oUXXlDt27c3rk9MTFSFCxdWQ4YMMSlfoUIF1aRJE+PryMhIZWdnp/bv329SLu17WrdunUksHh4exus1zeDBg5Wnp2em55GRtOumSpUqJp/DjBkzFKDatm1rUn7o0KEKMP4cHT9+XAFq4MCBJuV+++03Baj33ntPKaXUzZs3lbOzs8nno5RSu3fvVoDJtWzJ51GQE6Xces2HGjVqRKlSpfj666/566+/2L9//2Nvu27ZsoVChQrRsWNHk/VptyDTbmFu3boVgNdee82kXNeuXU1e379/n82bN9O+fXtcXV1JSUkxLmFhYdy/f599+/Zl6bxeeeUVk9d79uzhxo0b9OjRw+Q4er2eli1bsn//fuPto9q1a7NgwQImTZrEvn370Ol0WYrBUlWrVsXR0ZHXX3+dhQsXcu7cObO20+v1JueUmpr61HH89/m0s7MzZcqU4cKFC+nKPvo5//rrr6SkpBAREWESk7OzM40aNTLe4j1z5gwnTpwwXiOPfvexsbGcPHnS4tjnzp1L9erVcXZ2xt7eHgcHBzZv3szx48eNZQoXLkyvXr1YsGCB8TvfsmULf//9N4MHDzaWW7t2LRUrVqRq1aom8bVo0SLD29VNmzalSJEiJutq167NrVu36NKlCz/99BPx8fEWnU9YWBh2dg9/9ZYvXx6A1q1bm5RLWx8TEwM8/Bl89PFA7dq1KV++vPFnde/evdy/fz/dz2r9+vUJDg42WWfp51FQSaLMhzQaDb169eLbb79l7ty5lClThoYNG2ZY9vr16xQtWhSNRmOy3s/PD3t7e65fv24sZ29vj7e3t0m5okWLpttfSkoKs2bNwsHBwWQJCwsDsPgXS5qAgACT11evXgWgY8eO6Y710UcfoZTixo0bgOH5Wo8ePZg3bx716tXDy8uLiIgIrly5kqVYzFWqVCk2bdqEn58fgwYNolSpUpQqVYqZM2dmut3EiRNNzqdUqVJPFcej3xuAk5NTho1THvc516pVK93nvGzZMuP3mVZu+PDh6coNHDgQsPy7nzZtGgMGDKBOnTr8+OOP7Nu3j/3799OyZct0sb/55pvcvn2bJUuWAPDZZ59RvHhxXn75ZZNzOXLkSLr4ChcujFIqXXyPfhYA3bt35+uvv+bChQu88sor+Pn5UadOHaKiosw6Jy8vL5PXjo6Oma6/f/8+gPFnMaOYihUrZvKzCul/NjNaZ+nnUVDZ2zoAkT169uzJ2LFjmTt3Lh9++OFjy3l7e/Pbb7+hlDJJlnFxcaSkpODj42Msl5KSwvXr101+6T6aaIoUKYJWq6V79+4mDSj+KzQ0NEvn9GgyT4tt1qxZ1K1bN8Nt/P39jWVnzJjBjBkziImJYc2aNYwaNYq4uDg2bNgAGGpZjzZgAsMv97RjZUXDhg1p2LAhqampHDhwgFmzZjF06FD8/f159dVXM9zm9ddfN+n76OTklOXjW+pxn/OKFSvS1UgyKjd69Gg6dOiQYZmyZctaFMu3335L48aNmTNnjsn627dvpytbunRpWrVqxeeff06rVq1Ys2YNEyZMQKvVmsTo4uKSriHTo+eQ5tHPIk2vXr3o1asXd+/eZceOHYwbN442bdpw6tSpTD+jp5H2cxcbG0vx4sVN3rt8+bLJzyqk/9lMWxcSEmJ8bennUVBJosynAgMDGTFiBCdOnKBHjx6PLffCCy+wfPlyVq9eTfv27Y3rFy1aZHwfoEmTJkydOpUlS5bw1ltvGct99913JvtzdXWlSZMmHD58mMqVKxv/Ks4ODRo0wNPTM93ttScpUaIEgwcPZvPmzezevdu4PiQkhCNHjpiUPXXqFCdPnnziL4y0RJZZ9wGtVkudOnUoV64cS5Ys4dChQ49NlMWKFaNYsWLmnpJZx8+qFi1aYG9vz9mzZ9Pdlv2vsmXL8swzz/Dnn38yefJkqxxbo9Gk+yPhyJEj7N27l6CgoHTlhwwZQvPmzenRowdarZZ+/fqZvN+mTRsmT56Mt7d3lv9g+69ChQrRqlUrkpOTadeuHceOHcu2RNm0aVPA8MdDrVq1jOv379/P8ePHGTNmDAB169bF2dmZJUuWmHxfe/bs4cKFCyaJ0tqfR34liTIfmzJlyhPLRERE8Pnnn9OjRw/Onz9PpUqV2LVrF5MnTyYsLIwXX3wRgObNm/P8888zcuRI7t69S82aNdm9ezeLFy9Ot8+ZM2fy3HPP0bBhQwYMGEBISAi3b9/mzJkz/Pzzz2zZssUq5+fm5sasWbPo0aMHN27coGPHjvj5+XHt2jX+/PNPrl27xpw5c0hISKBJkyZ07dqVcuXKUbhwYfbv38+GDRtMaj7du3enW7duDBw4kFdeeYULFy4wdepUfH19nxhLpUqVjOfeo0cPHBwcKFu2LEuWLGHLli20bt2aEiVKcP/+feNf72mfrTU87vjWGGAiJCSEiRMnMmbMGM6dO0fLli0pUqQIV69e5ffff6dQoUJMmDABgC+++IJWrVrRokULevbsSWBgIDdu3OD48eMcOnSIH374waJjt2nThg8++IBx48bRqFEjTp48ycSJEwkNDc2wH22zZs2oUKECW7dupVu3bvj5+Zm8P3ToUH788Ueef/553n77bSpXroxerycmJoaNGzfyzjvvUKdOnUxj6tevHy4uLjRo0ICAgACuXLlCZGQkHh4eJgnM2sqWLcvrr7/OrFmzsLOzo1WrVpw/f57333+foKAg3n77bcBwV2f48OFMmjSJvn370qlTJy5evMj48ePT3Xq1xudRINi2LZGwlv+2es3Mo61elVLq+vXrqn///iogIEDZ29ur4OBgNXr0aHX//n2Tcrdu3VK9e/dWnp6eytXVVTVr1kydOHEiw5aK0dHRqnfv3iowMFA5ODgoX19fVb9+fTVp0iSTMljQ6vWHH37I8P3t27er1q1bKy8vL+Xg4KACAwNV69atjeXv37+v+vfvrypXrqzc3d2Vi4uLKlu2rBo3bpy6e/eucT96vV5NnTpVlSxZUjk7O6uaNWuqLVu2mNXqVSmlRo8erYoVK6bs7OwUoLZu3ar27t2r2rdvr4KDg5WTk5Py9vZWjRo1UmvWrMn0nLMio+MrZWit2Lp163TlHz2vJ11Dq1evVk2aNFHu7u7KyclJBQcHq44dO5p0w1FKqT///FN17txZ+fn5KQcHB1W0aFHVtGlTNXfu3Ceew6PX0oMHD9Tw4cNVYGCgcnZ2VtWrV1erV69WPXr0UMHBwRnuY/z48QpQ+/bty/D9O3fuqP/973+qbNmyytHR0diV6e2331ZXrlwxiWXQoEHptl+4cKFq0qSJ8vf3V46OjqpYsWKqc+fO6siRI5meW9p18/HHH5usf9z1ndH3kZqaqj766CNVpkwZ5eDgoHx8fFS3bt3UxYsXTbbV6/UqMjJSBQUFKUdHR1W5cmX1888/p/vOLfk8CnKrV41SMuyGECL/qFmzJhqNhv3799s6FJFPyK1XIUSel5iYyNGjR1m7di0HDx5k1apVtg5J5COSKIUQed6hQ4do0qQJ3t7ejBs3jnbt2tk6JJGPyK1XIYQQIhMy4IAQQgiRCUmUQgghRCYkUQohhBCZKHCNefR6PZcvX6Zw4cKPHZ5KCCFE/qeU4vbt2xQrVsxkoPqMCtrM9u3bjdPvAGrVqlVP3Gbbtm2qevXqysnJSYWGhqo5c+ZYdMyLFy8qQBZZZJFFFlkUkG7AhkfZtEZ59+5dqlSpQq9evTIdQzJNdHQ0YWFh9OvXj2+//Zbdu3czcOBAfH19zdoeMA7pdfHiRdzd3Z8qfpG76HQ6Nm7cSPPmzXFwcLB1OKIAkWsvb0pMTCQoKOiJQz3aNFG2atWKVq1amV1+7ty5lChRghkzZgCG+doOHDjAJ598YnaiTLvd6u7uLokyn9HpdLi6uuLu7i6/rESOkmsvb3vSY7g89Yxy7969NG/e3GRdixYtmD9/PjqdLsML9MGDByZTJyUmJgKGCzunJu8VOSPt+5TvVeQ0ufbyJnO/rzyVKK9cuWKcXzCNv78/KSkpxMfHZzihaWRkpHFmg//auHEjrq6u2RarsB1zJ9AVwtrk2stbkpKSzCqXpxIlpK8iq38HFnpc1Xn06NEMGzbM+DrtnnTz5s3l1ms+o9PpiIqKolmzZnL7S+QoufbyprQ7jE+SpxJl0aJF083aHRcXh729vXFW70c5OTllODu8g4ODXND5lHy3wlZy4trT6/UkJydn6zHyCwcHB7RababvmyNPJcp69erx888/m6zbuHEjNWvWlF+MQoh8Lzk5mejoaPR6va1DyTM8PT0pWrToU/Wbt2mivHPnDmfOnDG+jo6O5o8//sDLy4sSJUowevRoLl26xKJFiwDo378/n332GcOGDaNfv37s3buX+fPns3Tp0hyPPTbhHtHxdwn1KUSAh0uOH18IUbAopYiNjUWr1RIUFJR5B3mBUoqkpCTi4uIAMmzDYi6bJsoDBw7QpEkT4+u0Z4k9evRgwYIFxMbGEhMTY3w/NDSUdevW8fbbb/P5559TrFgx/u///s/sriHWsmx/DKNX/oVegZ0GIjtUIrxWiRyNQQhRsKSkpJCUlESxYsWkIaKZXFwMlZi4uDj8/PwyvQ2bGZsmysaNGxsb42RkwYIF6dY1atSIQ4cOZWNUmYtNuGdMkgB6Be+tPMrzZXylZimEyDapqakAODo62jiSvCXtjwqdTpflRCl1dwtFx981Jsk0qUpxPt68ZsZCCPE0ZIxqy1jj85JEaaFQn0LYPfK5azUaQnzkVogQQuRHearVa24Q4OFCZIdKPFjzDreUK3HKixb1qhFw9yRoA8HVG/LAQ3ZpjCSEyA22bdtGkyZNuHnzJp6enrYOJ0OSKLMgvKov/LLx4YoD/y4AWkcoHADugeBeDDwCwb04ePxncSkCNrx9Io2RhBDCfJIos0KfCk3+B7cvw+0rkHgZbsfCnThITYZbFwzL4zgUAs8g8Agy/OtZ4t8l2PBvId9sS6TSGEkIISwjiTIrnNyg0Yj061N1DxNn4j+QcAkSL0HCPw//vXsNdHfh2gnDkhEHV0PSLBIMRUKgSCh4hRr+LRIM9ulHGjJXZo2RJFEKUTDk9KOXBw8eMGLECL7//nsSExOpWbMm06dPp1atWsYyu3fv5r333uPkyZNUqVKFefPmUalSJQAuXLjA4MGD2bVrF8nJyYSEhPDxxx8TFhaW7bGDJErr0jr8W0MMAupkXEZ3z5BAE2Lg1kW4FQMJ//57K8aQZHVJcO24YUlHY6iJeoWCV0nwLgVepcC7tCGp2mfedDytMdJ/k6U0RhKi4LDFo5eRI0fy448/snDhQoKDg5k6dSotWrQwGXBmxIgRzJw5k6JFi/Lee+/Rtm1bTp06hYODA4MGDSI5OZkdO3ZQqFAh/v77b9zc3LI15v+SRJnTHFzAp7RhyUjKA0PN82Y03DwPNy8Y/n/jvOHf5DuGJJsQA9HbTbfVaA01Tu9nwCdtKWNYCvkADxsjvbfyKKlKodVomNyhotQmhSgAbPHo5e7du8yZM4cFCxYY5x/+6quviIqKYv78+cZa5bhx42jWrBkACxcupHjx4qxatYrOnTsTExPDK6+8YqxhlixZMltifRxJlLmNvZOhluhdKv17SsHdeLhxDm6cNfx7/Qxc//f/yXf+fe8cnP7VdFsXL/AtB75lCPctx4tdShKtKUFgUDABnlKbFKIgsMWjl7Nnz6LT6WjQoIFxnYODA7Vr1+b48ePGRFmvXj3j+15eXpQtW5bjxw131d566y0GDBjAxo0befHFF3nllVeoXLlytsSbEUmUeYlGA26+hqXEI7d2lTI8H71+GuJPGxJo/CnDcusi3LsBMXsMC+D974KzJ/hVAP8K//77rOFfZ5mCTIj8xhaPXh43FaJS6omDAaS937dvX1q0aMEvv/zCxo0biYyM5NNPP+XNN9/MnqAfIYkyv9BowD3AsIQ+b/pecpIhcV47CfEnIe64oSHRjXNw/5ZJAjXyDAb/ilC0IhStZFg8g23arUUI8XRs8eildOnSODo6smvXLrp27QoYhpM7cOAAQ4cONZbbt28fJUoYnpXevHmTU6dOUa5cOeP7QUFB9O/fn/79+zN69Gi++uorSZTCihxdIaCyYfkv3f2HiTPub7j6N1w9Zuj2ktbF5eQvD8s7efy7nyoPF+/SYJe18ROFEDkvvFYJni/jy/n4JEJ8XLO9fUKhQoUYMGAAI0aMMM4MNXXqVJKSkujTpw9//vknABMnTsTb2xt/f3/GjBmDj48P7dq1A2Do0KG0atWKMmXKcPPmTbZs2UL58uWzNe7/kkRZkDk4P0x4/5V0w5Awrx6FK0fhyhFDMn2QAOd3GpY0jm6G7YtVMyyBNQytb6XmKUSuFeDhkqMN+KZMmYJer6d79+7cvn2bmjVr8uuvv1KkSBGTMkOGDOH06dNUqVKFNWvWGAeAT01NZdCgQfzzzz+4u7vTsmVLpk+fnmPxa1Rm03fkQ4mJiXh4eJCQkIC7uzyHM1tKsuF27ZUjEPunYbnyl6Ery6NcvCCwOhSvBcVrGpKnS5H05axMp9Oxbt06wsLCZCJvkaNy4tq7f/8+0dHRhIaG4uzsnC3HyI8y+9zMzQdSoxTmsXd8ePu2WjfDutQUQ2Oh2D/g0iG4fMiQPO/dgDObDEsanzIQVBuC6hgW72fyxJi4QgghiVJkndbe0FrWvwJUNTykJ+WB4ZbtpUPwz37DcuPcwxa4h781lHMpAkF1oURdKFHPcNv2CYMlCCGELUiiFNZl72S41RpYA2r3A2DV7j/5Zd0aqmlOUcPuNDXto7G/dxNOrTcsAPYuEFQLghtAyHMQWNPwDFUIIWxMEqXIVrEJ93hn7T/oVXU2UR0AZ10qu3p643P9MMTsNSxJ1yF6h2EBsHc2POMMaQglGxkSr1aeOwohcp4kSpGtMhoJ5L7Sctq+HD71G0D9wYbBEq6dhAu74fwuw793rj5sYbttsmHGleD6UKoJlGwCfuWlZa0QIkdIohTZyqyRQDQa8CtnWGr1MSTO+NNwfgdE7zTUMu/dgDNRhgXAraghaZZ6wfDvv2PZCiGEtUmiFNkqSyOBaDTgW8aw1OoLer2hgVD0dji7FS7sgTtX4M+lhgUNFKuKXcmmFLnrZpgvFLlNK4SwDkmUIts99UggdnYPu6bUf9MwotDF3+DsZjizBa7+BZcPo718mOcBNfMzKP0ilGkJpV8AZ49sOS8hRMEgiVLkCKuOBOLgbGjgU7IRNJsIibFwdjP6UxtJPRWFQ9J1OLLMsNjZG55tlm0N5cLAM3vn3RNC5D+SKEXe5x4A1bqRWjGc9b+sIayiN/bnNsOpDYa+m2mtaTe8axjcvVwbw+L/rDQIEiKXaNy4MVWrVmXGjBm2DiUdixJlQkICq1atYufOnZw/f56kpCR8fX2pVq0aLVq0oH79+tkVpxBmURp7VHADKN0Ymn9gmKvz5Ho4uc7QDeXKX4ZlWyR4lYTyL0H5lw1D7knSFMLqHk2A27Zto0mTJty8eRNPT09juZUrV+baoSfNGkMsNjaWfv36ERAQwMSJE7l79y5Vq1blhRdeoHjx4mzdupVmzZpRoUIFli1blt0xC2E+71KGLii91sHwM/DybCjTCrROhhGDds+EeU1hZmXY+D+4dNDQ6lYIkaO8vLwoXLiwrcPIkFk1yipVqhAREcHvv/9OxYoVMyxz7949Vq9ezbRp07h48SLDhw+3aqBCPLVC3lDtNWJLduBCbBxlE/dS5MIGOLURbsXAnlmGpUgIVHzFsPg/a+uohcizevbsyfbt29m+fTszZ840eS9t5pAePXqwYMGCdDXPkJAQ+vbty6lTp1i5ciXe3t783//9H/Xr16dv375s3ryZ0NBQvvnmG2rWrJmt52FWojx27Bi+vr6ZlnFxcaFLly506dKFa9euWSU4Iaxt2f4YRq/8C70CO01hIjtMJPzl2Yb+mcdWG55r3jwPOz81LL7loXInqNRJGgKJ3EWpjGfvyQkOrmY9qpg5cyanTp2iYsWKTJw4kdTUVPbt20fHjh05efIk7u7uuLg8vpHf9OnTmTx5Mu+//z7Tp0+ne/fuNGjQgN69e/Pxxx/z7rvvEhERwbFjx9Bk46MTsxLlk5Lk05YXIifEJtwzJkkwDILw3sqjPF+mCQEVXoYKL0PyXUOyPLoSTm+Ea8dh80TDUqIeVHkVnm0vXU6E7emSYHIx2xz7vcvgWOiJxTw8PHB0dMTV1ZWiRYsC4O3tDYCfn5/JM8qMhIWF8cYbbwAwduxY5syZQ61atejUqRMA7777LvXq1ePq1avG/WcHixrzKKXYtGkTe/bs4cqVK2g0Gvz9/WnQoAEvvPBCtmZ0IZ5WRsPppSrF+fikh11XHAs9vO167xYc/9nQzeT8rofj0q5/F8qGGWZMKdUU7LQ5fi5CFASVK1c2/t/f3x+ASpUqpVsXFxeXOxLlpUuXaNOmDX/99RcVK1bE398fpRR79uzhgw8+MM5IHRgYmG3BCvE0zBpO779cPKF6d8OScAn++sEwEtC1E3BspWEpXMyQMKt2NTQcEiKnOLgaana2OnZOHOY/rWDTKmIZrdPr9dkah9mJcuDAgXh5eXHx4kUCAgJM3ouNjaVbt24MGjSI1atXWztGIawiS8PppfEIhOeGQoMhhomq/1hqSJy3L8POTwxLSEOo3sPQ5USmCBPZTaMx6/anrTk6OpKammryGjBZl9uZnSg3b97M7t270yVJgICAAD755BMaNmxo1eCEeJzYhHtEx98l1KeQRSP+PPVwehqNYZLpYtUM/TRProfDi+HM5oeznbgUgcqvQs1e4FvWwjMTIn8JCQnht99+4/z587i5uREcHIxGo2Ht2rWEhYXh4uKCm5ubrcPMlFn9KMHQqvXGjRuPff/mzZuZtl4SwlqW7Y+hwZQtdP3qNxpM2cKy/TEWbR/g4UK9Ut5PP6SevRM82w66/QhvH4XGo8G9ONy7Cb/Ngc9rw4I2hoZBKclPdywh8qjhw4ej1WqpUKECvr6+6HQ6JkyYwKhRo/D392fw4MG2DvGJNEqZ17v6zTff5KeffmLatGk0a9YMDw9Dq7+EhASioqJ45513aNeuXbq+MrlNYmIiHh4eJCQk4O7ubutwhIViE+7RYMqWdM8Zd41qgo+rPevWrSMsLMx2I3zoU+HsFjjwtaH1rPr32YmbP9ToZahlFs6+RgfCNnQ6XbZfe/fv3yc6OprQ0FCcneXWvrky+9zMzQdm33r99NNPSUlJ4bXXXiMlJcV4nzk5ORl7e3v69OnDxx9/nMVTEcI8mbVc9SmRC/7wsdPCM80My62LcGghHFpkmIh6+xTDs8wKL0OdARBUy9bRCiHMYHaidHR0ZM6cOXz00UccOHCAq1evAlC0aFFq1KghtTORIyxuuWpLnkHQ9H/w/Eg4vgZ+/wou7oOjPxqWwJpQd4AhcWpz5xiXQogszB7i7u5O06ZNsyMWIZ4os5arOp3O1uFlzN4RKnU0LLF/wm9fGFrMXjoAP/aBqLFQpz/U6CEDGQiRC5mdKP/v//7PrHJvvfVWloMRwhxP3XLVlgKqQLvZ8OJ4w3PM/fMh8RJEvQ/bpxqSZd0B4FHc1pEKIf5ldqKcPn26yeu0/pT29g93odFoLE6Us2fP5uOPPyY2NpZnn32WGTNmZNrNZMmSJUydOpXTp0/j4eFBy5Yt+eSTT4zDIomCwaoTQduCmx80HgUNhhpql3s/MwxksPcz+G0uVA439NmU7iXiEWa2vxT/ssbnZXb3kOjoaJPFxcWF7du3m6w7d+6cRQdftmwZQ4cOZcyYMRw+fJiGDRvSqlUrYmIybu6/a9cuIiIi6NOnD8eOHeOHH35g//799O3b16LjCpFrODgbRv4ZsBdeW2EYtECfAn8sMXQv+f41w9RfosDTag1DJSYnS1cjSyQlGQaOf5rWyBY/o7SmadOm0adPH2OimzFjBr/++itz5swhMjIyXfl9+/YREhJirLWGhobyxhtvMHXq1ByNWwirs7N72Fr2nwOwazqc+AVOrDUspZpCw+EQ0sDWkQobsbe3x9XVlWvXruHg4ICdndn1nAJJKUVSUhJxcXF4enoa/9DICpslyuTkZA4ePMioUaNM1jdv3pw9e/ZkuE39+vUZM2YM69ato1WrVsTFxbFixQpat2792OM8ePCABw8eGF8nJiYChn5Pubbxh8iStO8zz3+v/lXglQUQfwrtnplojq5Ac3YLnN2CvkQ99A1HooKfM2uaI5Ezcura8/X1JSYmhvPnz2frcfITd3d3vL29M/xuzP2+bJYo4+PjSU1NNY7+nsbf358rV65kuE39+vVZsmQJ4eHh3L9/n5SUFNq2bcusWbMee5zIyEgmTJiQbv3GjRtxdc2FXQrEU4uKirJ1CNZj3xrX8rUpHfcLJa7vQBuzF7sl7bnoXJa//dtzv0h5SZi5SE5de1qtVmZrMkNqamqmzyjTbss+idmJMq0mlkaj0XDnzp106y3tT/nol62UeuwF8Pfff/PWW28xduxYWrRoQWxsLCNGjKB///7Mnz8/w21Gjx7NsGHDTM4jKCiI5s2bS9/PfEan0xEVFUWzZs1sNzJPtumBPvEyZ36aTMiFFQTdP0nQhSlcTayBd5vxqBL1bB1ggZa/r73869H89ThmJ0pPT0+TBKaUolq1aiavNRqN2SPC+/j4oNVq09Ue4+Li0tUy00RGRtKgQQNGjBgBGOYqK1SoEA0bNmTSpEkZDtju5OSEk5NTuvUODg5yQedT+fW7jbX3o9Xptviq5xhgv4Yu2i343zwIi18yPMNs+j8IrGHrMAu0/Hrt5VfmfldmJ8qtW7dmOZiMODo6UqNGDaKiomjfvr1xfVRUFC+//HKG2yQlJZl0R4GHLcGkybTI79KG77uKF+NTevJFyksMtl9NF4ft2P37DJPybaHp++BbxtbhCpFvmJ0oHzx4QJMmTaz619KwYcPo3r07NWvWpF69enz55ZfExMTQv39/wHDb9NKlSyxatAiAl156iX79+jFnzhzjrdehQ4dSu3ZtihUrZrW4hMiNHh2+LxZvxqb2pVm/yfgdmgF/fm8YKu/EWsNE0o3fM8yjKYR4Kma3L+7fvz++vr6Eh4fz3XffcevWrac+eHh4ODNmzGDixIlUrVqVHTt2sG7dOoKDgwHDhND/7VPZs2dPpk2bxmeffUbFihXp1KkTZcuWZeXKlU8dixC5Xdrwfdp/H4GkDd/nF1wO2s+FAXugbGvDjCWHv4VZ1WHTeLifYNvAhcjjzJ5mC+DIkSOsWbOGNWvWcOTIERo0aMDLL79M27ZtCQkJycYwrUem2cq/cmKqo9wgNuFe5sP3XfwdosZBzL/drFy8oNG7UKuPDL6eTQrKtZffmJsPLOqxWrlyZf73v//x+++/c+7cOTp16sSGDRsoX748VapUYezYsRw4cOCpgxdCPN4TJ54Oqg291sGrS8GnDNy7ARvehdl14cQ6kOf5Qlgky0M7FCtWjP79+7Nu3TquXbvG+++/z/nz52nZsiWTJ0+2ZoxCCEtpNFAuzDA0XpvpUMgXrp+B77vAwpfgylFbRyhEnmGVAQfc3Nzo2LEjHTt2RK/Xc/36dWvsVgjxtLT2ULM3VOwIu6bB3tlwfid80RBq9IQm/4NCMqGAEJnJUqLcvHkzmzdvJi4uDr1eb1yv0WiYP38+vr6+VgtQCGEFzu6Gqb1q9IJN4+DYKsM0X0d/hCZjoGYfQ1IVQqRj8a3XCRMm0Lx5czZv3kx8fDw3b940Ljdu3MiOGIUQ1lIkGDotgJ6/gH8lQ4vY9SPhy0ZwIeMxloUo6Cz+E3Lu3LksWLCA7t27Z0c8QoicEPIcvLEdDi6ALR/A1aPwTSvDPJjNPoDCGY+OJURBZHGNMjk5mfr162dHLEKInGSnNXQZGXzQ8LwSDRxZBp/Vgt+/Ar1hOMrYhHvsORtPbMI9m4YrhK1YnCj79u3Ld999lx2xCCFsoZA3vDQT+m2GgKrwIAHWDYd5L7Bx03oaTNlC169+o8GULSzbn/Gk6kLkZ2bdev3v7Bt6vZ4vv/ySTZs2Ubly5XSda6dNm2bdCIUQOSOwBvTbYmjks3kiXD7MC5e68J62JdNSOpGknHlv5VGeL+P7+D6cQuRDZiXKw4cPm7yuWrUqAEePmvbFkvnRhMjj7LRQux+Uf4lrP76D7/m19LVfT0vtfv6n68U2fTXOxydJohQFilmJ0tozhwghcrnCRUlpP49eU6fxgf3XFNfEs8DxY9ak1qeka2VA+l6KgiPLI/MIIfK3AA8XWrbrTsvkj/kypTWpSkNb7R78FzeCv1bIUHiiwDArUfbv35+LFy+atcNly5axZMmSpwpKCJE7hNcqQdSoVlTqNYubXdeD37OQdB1+7ANLX4XEWFuHKES2M+vWq6+vLxUrVqR+/fq0bduWmjVrUqxYMZydnbl58yZ///03u3bt4vvvvycwMJAvv/wyu+MWQuSQAA+Xf59J1oNS22D3TNgxFU5tgNl1oOVHUOVVw/iyQuRDZtUoP/jgA06fPs3zzz/P3LlzqVu3LiVKlMDPz4+yZcsSERHBuXPnmDdvHnv37qVSpUrZHbcQwhbsHaHRCHhjBxSrZhjZZ3V/qV2KfM3skXn8/PwYPXo0o0eP5tatW1y4cIF79+7h4+NDqVKlpMWrEAWJX3noswl2z4BtU/6tXdaFNtOg4iu2jk4Iq8rSKMienp54enpaORQhRJ6itYfnh0PZMEOtMvZPWNHbMOdl2Mfg6mXrCIWwCmn1KoR4Ov4VoO9maPQuaLRwdAXMqQ9nt9g6MiGsQhKlEOLpaR2gyXvQJwq8n4HbsbC4PWwYDbr7to5OiKciiVIIYT3Faxga+tTqa3i9bzZ81QSuHM18OyFyMUmUQgjrcnSF1p9C1x+gkC/E/Q1fNYXfvgClZDYSkefIlOZCiOxRpjkM2AtrBhtaxa4fyaWDv/DyxS7EK3fsNBDZoRLhtUrYOlIhMmVWoqxWrZrZ3T8OHTr0VAEJIfIRN1/o8j38/hVq4/8IjNvOL45/MlQ3iL36Z2U2EpEnmJUo27VrZ/z//fv3mT17NhUqVKBevXoA7Nu3j2PHjjFw4MBsCVIIkYdpNFDndf7UVqDQmtd5xu4SSxwmMyu1Hf+X0kFmIxG5nlmJcty4ccb/9+3bl7feeosPPvggXRlzx4MVQhQ8/s/U4AXdJN7XLqKL/VaG2K+irt1xQpyXIrORiNzM4sY8P/zwAxEREenWd+vWjR9//NEqQQkh8p8ADxfGdajJ/1Jf563kQdxRztSxO4H/khfhzGZbhyfEY1mcKF1cXNi1a1e69bt27cLZ2dkqQQkh8qfwWiXYNaoJXfq8Q1KvLVC0kmE2km9fga2TQZ9q6xCFSMfiVq9Dhw5lwIABHDx4kLp16wKGZ5Rff/01Y8eOtXqAQoj85eFsJN6G8WI3jIKD38D2j+Dib9BhnqERkBC5hMWJctSoUZQsWZKZM2fy3XffAVC+fHkWLFhA586drR6gECIfc3CGl2ZAiXqwdiic2wZfPA+dF0JQbRsHJ4RBlvpRdu7cWZKiEMJ6qoRDQBVY3h3iT8E3YdBiMtTuJ/NcCpuTkXmEELmDXznotwUqtAO9DtaPgJX9IDnJpJiM7CNymsU1ytTUVKZPn87y5cuJiYkhOTnZ5P0bN25YLTghRAHjVBg6LTCMEbvxffjrB4g7Aa9+C0VCWLY/htEr/0KvkJF9RI6xuEY5YcIEpk2bRufOnUlISGDYsGF06NABOzs7xo8fnw0hCiEKFI0G6g2CHj8bxoq9+hd82ZjrRzYYkySAXsF7K49KzVJkO4sT5ZIlS/jqq68YPnw49vb2dOnShXnz5jF27Fj27duXHTEKIQqikAbw+nYoVh3u3cRrVRf62f0MKGORVKU4H5/0+H0IYQUWJ8orV65QqVIlANzc3EhISACgTZs2/PLLL9aNTghRsHkEQq/1UK0bGqVntMNSpjvMxgnDIx+tRkOIj6uNgxT5ncWJsnjx4sTGxgJQunRpNm7cCMD+/ftxcnKybnRCCOHgDG0/g7BP0Gu0tNfu5gfHCQRqbjC5Q0UZJ1ZkO4sTZfv27dm82TDc1JAhQ3j//fd55plniIiIoHfv3lYPUAgh0Gigdj/sIn5C7+xFZbtodhSZQHjRK7aOTBQAFrd6nTJlivH/HTt2pHjx4uzZs4fSpUvTtm1bqwYnhBAmQhti98Y2+L4r2qtHYUEbaDvL0A9TiGzy1BM3161b1ziUnRBCZLsiwdD7V1j1BpxYC6teh2vHoelYsJOu4cL6snRVLV68mAYNGlCsWDEuXLgAwIwZM/jpp5+sGpwQQmTIyQ06L4aG7xhe75oOy7pB8l3bxiXyJYsT5Zw5cxg2bBhhYWHcunWL1FTDaP+enp7MmDHD4gBmz55NaGgozs7O1KhRg507d2Za/sGDB4wZM4bg4GCcnJwoVaoUX3/9tcXHFULkcXZ28MJY6PAVaJ3g5C/wdUtIvGzryEQ+Y3GinDVrFl999RVjxoxBq9Ua19esWZO//vrLon0tW7aMoUOHMmbMGA4fPkzDhg1p1aoVMTExj92mc+fObN68mfnz53Py5EmWLl1KuXLlLD0NIUR+UbmzYXACVx+4cgS+agqX/zC+LUPeiadl8TPK6OhoqlWrlm69k5MTd+9adttj2rRp9OnTh759+wKG27e//vorc+bMITIyMl35DRs2sH37ds6dO4eXlxcAISEhlp6CECK/KVEH+m2G78Lh2gn4phW8Mp9ltyvKkHfiqVmcKENDQ/njjz8IDg42Wb9+/XoqVKhg9n6Sk5M5ePAgo0aNMlnfvHlz9uzZk+E2a9asoWbNmkydOpXFixdTqFAh2rZtywcffICLS8Z9qR48eMCDBw+MrxMTEwHQ6XTodDqz4xW5X9r3Kd9rAeUWCBHr0K7qg925rajvu3IipTt61RIwDHk3euVf1AstQoCHdSeZl2svbzL3+7I4UY4YMYJBgwZx//59lFL8/vvvLF26lMjISObNm2f2fuLj40lNTcXf399kvb+/P1euZNw36ty5c+zatQtnZ2dWrVpFfHw8AwcO5MaNG499ThkZGcmECRPSrd+4cSOurjKiR34UFRVl6xCEDWncu1PZWxFyfRvj7BcRRByTUrqhxw69guXrtvKMh3ryjrJArr28JSnJvOEPLU6UvXr1IiUlhZEjR5KUlETXrl0JDAxk5syZvPrqqxYHqnlkrjmlVLp1afR6PRqNhiVLluDh4QEYbt927NiRzz//PMNa5ejRoxk2bJjxdWJiIkFBQTRv3hx3d3eL4xW5l06nIyoqimbNmuHg4GDrcIQtqZe4uXUaRfZG0tt+A8U113hLN5hkjROdw5pkS41Srr28J+0O45NkqR9lv3796NevH/Hx8ej1evz8/Czeh4+PD1qtNl3tMS4uLl0tM01AQACBgYHGJAlQvnx5lFL8888/PPPMM+m2cXJyynBoPQcHB7mg8yn5bgVAkRaj2HPfixqH3qO59iBLNR9yofl8SvgUzrZjyrWXt5j7XT1V71wfH58sJUkAR0dHatSoke5WRVRUFPXr189wmwYNGnD58mXu3LljXHfq1Cns7OwoXrx4luIQQuRf9V9+ndvhK0hx9KCa3RnaHeoF18/aOiyRx1icKK9evUr37t0pVqwY9vb2aLVak8USw4YNY968eXz99dccP36ct99+m5iYGPr37w8YbptGREQYy3ft2hVvb2969erF33//zY4dOxgxYgS9e/d+bGMeIUTB5lOhMfb9NoFHCbhxDuY3g38O2joskYdYfOu1Z8+exMTE8P777xMQEPDY54nmCA8P5/r160ycOJHY2FgqVqzIunXrjC1qY2NjTfpUurm5ERUVxZtvvknNmjXx9vamc+fOTJo0KcsxCCEKAN8y0HcTfNcJYv+EhW2g8yJ4ppmtIxN5gEYpZVHzr8KFC7Nz506qVq2aTSFlr8TERDw8PEhISJDGPPmMTqdj3bp1hIWFyXMikbEHd2B5dzi7BTRaePkzqNr1qXcr117eZG4+sPjWa1BQEBbmViGEyB2c3KDLMqgcDioVVg+AnZ/Cv7/TZBQfkRGLb73OmDGDUaNG8cUXX8ioOEKIvMfeEdrNhcJFYfdM2DwR7lxjmXd/Rq86JqP4iHTMSpRFihQxeRZ59+5dSpUqhaura7rbDDdu3LBuhEIIYW12dtBsIrgVhV9Hw29zcE79E63qjx579AreW3mU58v4EuAhDQULOrMSZVZmBRFCiFyv3kAo5IN+1QBe1u7BkzsM0A0lCWdSleJ8fJIkSmFeouzRo0d2xyGEELZRuTO3VCGcV/aikfYISzST6Zk8kjuawoT4yDCX4ikHHBBCiPzAq0pr9jSYzy1ViGp2Z/jBcSLTwvykNikASZRCCAHAi81fQhfxC8ku/pSx+4eXD8ooPsJAEqUQQvzLt1Q1HF/fCEVC4VYMfN0Srh6zdVjCxiRRCiHEfxUJgd6/gn9FuBsH34TBPwdsHZWwIUmUQgjxqML+0HMtFK8N92/BwrYQvcPWUQkbMavVa4cOHcze4cqVK7McjBBC5BouRaD7Klj2GpzbBt92hM4LoWwrW0cmcphZNUoPDw/j4u7uzubNmzlw4OGtiIMHD7J582aTeSKFECLPSxvyrlwbSH0Ay7rB0R9tHZXIYWbVKL/55hvj/9999106d+7M3LlzjdNqpaamMnDgQBlkXAiR/zg4Q6eF8NNAOLIMVvSB5CSo3t3WkYkcYvEzyq+//prhw4ebzD2p1WoZNmwYX3/9tVWDE0KIXEFrbxgftkYvQMGawbBvrvHt2IT7nE7QEJtw33YximxjcaJMSUnh+PHj6dYfP34cvV5vlaCEECLXsbODNtOh7iDD6w3vws5pLNsfQ+NPd/DZ31oaf7qDZftjMt+PyHMsnj2kV69e9O7dmzNnzlC3bl0A9u3bx5QpU+jVq5fVAxRCiFxDo4EWH4JjIdgxFTZP4ErKX+jVK4BGBlPPpyxOlJ988glFixZl+vTpxMbGAhAQEMDIkSN55513rB6gEELkKhoNNB1jeHa5eSJD7FfihI4pKa8CGhlMPR+yOFHa2dkxcuRIRo4cSWJiIoA04hFCFDwN3yEhxR6P7WPpb/8zTiQzISUCrcZOBlPPZyxOlGB4Trlt2zbOnj1L165dAbh8+TLu7u64ublZNUAhhMitPJoM4cB1HTWPfkAv+19x1KRg/9I0qU3mMxYnygsXLtCyZUtiYmJ48OABzZo1o3DhwkydOpX79+8zd+7cJ+9ECCHyiZodh3PdzxOvLcN5TbsZYj8G/f+BnfbJG4s8weJWr0OGDKFmzZrcvHkTF5eHfzW1b9+ezZs3WzU4IYTIC9zr9eBQ8OsojR0c/hZWD4DUFFuHJazE4hrlrl272L17N46Ojibrg4ODuXTpktUCE0KIvOQfrwZUqV4L+9VvGAYm0KdA+y8NfTBFnmZxjVKv15Oamppu/T///EPhwoWtEpQQQuRFqkI7w3iwdvaGoe5+7AOpOluHJZ6SxYmyWbNmzJgxw/hao9Fw584dxo0bR1hYmDVjE0KIvKf8S9B5Mdg5wN+rYUUvSEm2dVTiKVicKKdPn8727dupUKEC9+/fp2vXroSEhHDp0iU++uij7IhRCCHylnJh8OoS0DrC8Z/hh57GZBmbcI89Z+OJTbhn2xiF2Sy+eV6sWDH++OMPli5dyqFDh9Dr9fTp04fXXnvNpHGPEEIUaGVawKtL4fuucPIX+KEHP5ScxLurT6BXYKeByA6VCK9VwtaRiifI0lNmFxcXevfuTe/eva0djxBC5B/PvAhd0pLlOor8fQV7NYRkHGS4uzzE4luvWq2WJk2acOPGDZP1V69eNZlRRAghBFD6BeiylFStEy9qDzHHYQaOGBr4pA13J3I3ixOlUooHDx5Qs2ZNjh49mu49IYQQjyjVlFvtvuWecuQF7WFjstRqNDLcXR5gcaLUaDT8+OOPvPTSS9SvX5+ffvrJ5D0hhBDpeVdqzr66s43Jcq7DDKa8XEZuu+YBWapRarVaZs6cySeffEJ4eDiTJk2S2qQQQjxBk1adSOr4HalaZ5pqD9Pp7HuQ8sDWYYkneKohI15//XXKlClDx44d2b59u7ViEkKIfMu7UjMotBy+C4fTv8Ky7hC+GOydbB2aeAyLa5TBwcEmjXYaN27Mvn37+Oeff6wamBBC5FslG0HXZWDvYkiWyyOkZpmLWZwoo6Oj8fb2NllXunRpDh8+zLlz56wWmBBC5GslG0HX78HeGU5tgOU9ZASfXMriRPk4zs7OBAcHW2t3QgiR/5VsDF3SkuV67n/Xjb2nLsuoPbmMWYnSy8uL+Ph4AIoUKYKXl9djFyGEEBYo1cTQz9LOEedzv5K4uDuNpmxk2f4YW0cm/mVWY57p06cbZwb574DoQgghnl6sTz3evT+Mrxw+pYX2ADOZxdCVb8moPbmEWYmyR48eGf5fCCHE04uOv8sOfWVe1w3jS4dPaaXdj55ZXIirTYCHC7EJ94iOv0uoTyFJnDZgVqJMTEw0e4fu7u5ZDkYIIQqiUJ9C2Glgu74Kb+je5guH6bTW/s69/SNYfmMCo1Yfl4HUbcisZ5Senp4UKVIk0yWtjKVmz55NaGgozs7O1KhRg507d5q13e7du7G3t6dq1aoWH1MIIXKTAA8XIjtUQqvRsE1fjUEpb5Oqscfl1Bqc1w5Co1IBjAOpS2OfnGVWjXLr1q3ZcvBly5YxdOhQZs+eTYMGDfjiiy9o1aoVf//9NyVKPP4vpoSEBCIiInjhhRe4evVqtsQmhBA5KbxWCZ4v48v5+CRCfJqija2OfnkEbdlDKhre0Q1Aj51xIHW5BZtzzEqUjRo1ypaDT5s2jT59+tC3b1/A0FDo119/Zc6cOURGRj52uzfeeIOuXbui1WpZvXp1tsQmhBA5LcDD5WEC9AjjVusvKLymH+21u9FjxwjdG2g0WhlIPYdleQi7pKQkYmJiSE427SBbuXJls7ZPTk7m4MGDjBo1ymR98+bN2bNnz2O3++abbzh79izffvstkyZNeuJxHjx4wIMHD0e8SHveqtPp0Ol0ZsUq8oa071O+V5HTsuvaK1y5Lfsu3qDe4ZG8ot2JHjt0YdPxcbWX69wKzP0MLU6U165do1evXqxfvz7D91NTU83aT3x8PKmpqfj7+5us9/f358qVKxluc/r0aUaNGsXOnTuxtzcv9MjISCZMmJBu/caNG3F1lb/K8qOoqChbhyAKqGy59jR+7Co+gOcvzaGTdjvnD7zNuqu9QGO18WIKrKQk8+YCtThRDh06lJs3b7Jv3z6aNGnCqlWruHr1KpMmTeLTTz+1ONBHp+ZSSmU4XVdqaipdu3ZlwoQJlClTxuz9jx49mmHDhhlfJyYmEhQURPPmzaWFbj6j0+mIioqiWbNmODg42DocUYBk/7UXhv5YZTQ/DSDk+naCgkPRt/wYZGrDp2Jujw6LE+WWLVv46aefqFWrFnZ2dgQHB9OsWTPc3d2JjIykdevWZu3Hx8cHrVabrvYYFxeXrpYJcPv2bQ4cOMDhw4cZPHgwAHq9HqUU9vb2bNy4kaZNm6bbzsnJCSen9KPyOzg4yC/TfEq+W2Er2XrtVX3VkBhXvYH20AK0WgcIk2T5NMz9riyuu9+9exc/Pz/AMLTdtWvXAKhUqRKHDh0yez+Ojo7UqFEj3a2KqKgo6tevn668u7s7f/31F3/88Ydx6d+/P2XLluWPP/6gTp06lp6KEELkLVXCod1sQAP7v4INo0DmAs52Ftcoy5Yty8mTJwkJCaFq1ap88cUXhISEMHfuXAICAiza17Bhw+jevTs1a9akXr16fPnll8TExNC/f3/AcNv00qVLLFq0CDs7OypWrGiyvZ+fH87OzunWCyFEvlW1K+hTYM2b8Ntc0GihxYdSs8xGWXpGGRsbC8C4ceNo0aIFS5YswdHRkQULFli0r/DwcK5fv87EiROJjY2lYsWKrFu3zjgLSWxsLDExMjCwEEKYqB4BSg8/D4F9n3NHp+dI+XcI9XWT/pXZQKPU09Xbk5KSOHHiBCVKlMDHx8dacWWbxMREPDw8SEhIkMY8+YxOp2PdunWEhYXJM0qRo2x27e2fD78YGivOTmnLJ6nhRHaoLEPcmcncfPDU7YtdXV2pXr16nkiSQgiRn8SW6co4XU8ABtqv4W3tD7y38i8Z4s7KLL71qpRixYoVbN26lbi4OPR6vcn7K1eutFpwQgghHi86/i4LU5ujQc94h0W8ab8aPXacj68rt2CtyOIa5ZAhQ+jevTvR0dG4ubnh4eFhsgghhMgZabOOLEhtyURddwCG2K+k4unZNo4sf7G4Rvntt9+ycuVKwsLCsiMeIYQQZkqbdeS9lUf5OrUV9ho979kvofC+T8DFCRqNtHWI+YLFidLDw4OSJUtmRyxCCCEs9OisIxwtC1FjYauhy0hslcEy6fNTsvjW6/jx45kwYQL37snDYiGEyA0CPFyoV8rbkAgbDIEXxxve2DKJbz8eQtevfqPBlC0s2y/d7bLC4hplp06dWLp0KX5+foSEhKRrCm3J6DxCCCGywXNvk3gvGffdkxlhvwy90jAntS3vrTzK82V8pWZpIYsTZc+ePTl48CDdunXD398/wwHMhRBC2NbRkn3Yve0MIxyW867D9yhgbmpbmfQ5CyxOlL/88gu//vorzz33XHbEI4QQwgpCfQrRTd8OO52edxxWMMrhezQajeE5prCIxc8og4KCZEQbIYTI5dJaxM7Wv8Knuo4AvGu/lIC/5to4srzH4kT56aefMnLkSM6fP58N4QghhLCW8Fol2DWqCfV7T+V2vX+7imwaD7um2zSuvMbiW6/dunUjKSmJUqVK4erqmq4xz40bN6wWnBBCiKcT4OFieCZZagw4Oxi6jWwaD0pPbOWB0nXEDBYnyhkzZmRDGEIIIbJdo5GABrZOgs0TWfLrCT5LaYedBiI7VJLB1B/DokSp0+nYtm0b77//vgw6IIQQeVGjESQmp+K+O5Lh9stBKT5LbS9dRzJh0TNKBwcHVq1alV2xCCGEyAFHS/blI92rAAx3+IEh2h9JVXrOxyfZOLLcyeLGPO3bt2f16tXZEIoQQoicEOpTiC/0bZms6wLA2w4/Mtx+BSHeUpvMiMXPKEuXLs0HH3zAnj17qFGjBoUKFTJ5/6233rJacEIIIazv4WDqGlJ1drzvsITB9qvgQAi8MA5kIBkTFifKefPm4enpycGDBzl48KDJexqNRhKlEELkAQ8HU69DwsXyeGz7n6HbSKoOmk+SZPkfFifK6Ojo7IhDCCFEDnvYdeRNcHWGdcNh72eQ8gBaTQU7i5/O5UtP9SkopVBKWSsWIYQQtlK7H7w0E9DA/q9g7RDQ620dVa6QpUS5aNEiKlWqhIuLCy4uLlSuXJnFixdbOzYhhBA5qUZPaDcHNHZwaBFx3/Yi9uZtW0dlcxYnymnTpjFgwADCwsJYvnw5y5Yto2XLlvTv35/p02VYJCGEyNOqdmFP1Y9IUXb4nVvNoWmvsPy3s7aOyqYsfkY5a9Ys5syZQ0REhHHdyy+/zLPPPsv48eN5++23rRqgEEKInBObcI9u+wJ5UTOEzxz+j9ba39i0ti+xpVcS4F3E1uHZhMU1ytjYWOrXr59uff369YmNjbVKUEIIIWwjOv4uegUb9bXopxvOfeXAi9pDuKzoBslJxCbcY8/ZeGIT7tk61BxjcaIsXbo0y5cvT7d+2bJlPPPMM1YJSgghhG2E+hTC7t+eIdv1Veipe5e7ygnP2F1cm9OallPW0vWr32gwZQvL9sfYNtgcYvGt1wkTJhAeHs6OHTto0KABGo2GXbt2sXnz5gwTqBBCiLzj4WAER0lViv3qWfY2mEfTg4PwvXmIbx0mEZE8ipvKvcCMD2txonzllVf47bffmD59OqtXr0YpRYUKFfj999+pVq1adsQohBAiBz0cjCCJEB9XAjxc+LOIB4Fru1HJ7jzLHT+gW/JoriovzscnSaLMSI0aNfj222+tHYsQQohcwjgYwb/8ytTmVd37LHKI5Bm7S6xwnECEbgwhPq42jDJnyLALQgghnijAw4V+7VsSnjyeaL0/QXbXWOc+iYB7+b/riNmJ0s7ODq1Wm+lib5+lCqoQQog8ILxWCZaP6sz1zmvQ+VTA5UE8fBMGMftsHVq2MjuzZTYP5Z49e5g1a5YMZyeEEPmc4ZZsOSi1Hr4Lh4v7YFE7brSZx4nCdQn1KZTvnlmanShffvnldOtOnDjB6NGj+fnnn3nttdf44IMPrBqcEEKIXMrFE7qvguURcCaKwqsiWK57gzXqOSI7VCK8VglbR2g1WXpGefnyZfr160flypVJSUnhjz/+YOHChZQokX8+GCGEEE/g6Eps2Hx+Sq2PgyaVGY6z6WX3C++tPJqvBiSwKFEmJCTw7rvvUrp0aY4dO8bmzZv5+eefqVixYnbFJ4QQIheLvqljqG4g81NaAfC+wxJGaJdw/tpdG0dmPWYnyqlTp1KyZEnWrl3L0qVL2bNnDw0bNszO2IQQQuRyoT6F0Gjs+CClG1N0rwLQ334tVQ68CynJNo7OOsx+Rjlq1ChcXFwoXbo0CxcuZOHChRmWW7lypdWCE0IIkbv9dySfualtuYEnUxy/wvXECh4sisOp6xJwdrd1mE/F7EQZERGBRqPJzliEEELkQf8dyefIpXL03uDB5w4zKBSzg5uzX6RI35/APcDWYWaZ2YlywYIF2RiGEEKIvCytS8hr8/ahV1UIT36fbxw/xjfxJKlfvcCNdt9ymhJ5svuIjMwjhBDCKtKm6AI4qkrSPnkCZ/UBaG9fwnlRGJ/Pn5cnZx2RRCmEEMIq/jtFF8A/yo+OyRP4XV+Owpp7LHCYyit22/Jc9xGbJ8rZs2cTGhqKs7MzNWrUYOfOnY8tu3LlSpo1a4avry/u7u7Uq1ePX3/9NQejFUII8ThpDXu0/7Zn0Wo0dGpYiW7Jo1n9b1/Ljx2+5B3tUs5fu2PjaM1n00S5bNkyhg4dypgxYzh8+DANGzakVatWxMRkXC3fsWMHzZo1Y926dRw8eJAmTZrw0ksvcfjw4RyOXAghREbCa5Vg16gmLO1Xl12jmtDruVBSNA4M1Q1iZkp7AAbar6HqviGQnDf6WmqUDQdorVOnDtWrV2fOnDnGdeXLl6ddu3ZERkaatY9nn32W8PBwxo4da1b5xMREPDw8SEhIwN09bzdZFqZ0Oh3r1q0jLCwMBwcHW4cjChC59jK3bH+McSLoDtpdfOz0FVq9DgKqQJfvwb2YTeIyNx/YbLqP5ORkDh48yKhRo0zWN2/enD179pi1D71ez+3bt/Hy8npsmQcPHvDgwQPj68TERMBwYet0uixELnKrtO9TvleR0+Tay1yHqgHUCy1CzI0kSng9j0psjVrRA03sn6gvG5PacREqsEaOx2Xu92WzRBkfH09qair+/v4m6/39/bly5YpZ+/j000+5e/cunTt3fmyZyMhIJkyYkG79xo0bcXXN/xOOFkRRUVG2DkEUUHLtPdl14DDgEjKauuem437nHzQL27DRpzfXfBvg6ZRzsSQlJZlVzuYTSD46iIFSyqyBDZYuXcr48eP56aef8PPze2y50aNHM2zYMOPrxMREgoKCaN68udx6zWd0Oh1RUVE0a9ZMbn+JHCXXXhY96MDFRb0IittG2LUvmBd7kaSwiXSsFZIjh0+7w/gkNkuUPj4+aLXadLXHuLi4dLXMRy1btow+ffrwww8/8OKLL2Za1snJCSen9H+iODg4yAWdT8l3K2xFrj3LxCa50PhiX97SFmGI/Sr62q9j14YYrpdcRtGixbP9+OZ+VzZr9ero6EiNGjXS3aqIioqifv36j91u6dKl9OzZk++++47WrVtnd5hCCCGySXT8XVKVHdNTOjEw+S2SlBPP2R2lyOJmcDn39GawafeQYcOGMW/ePL7++muOHz/O22+/TUxMDP379wcMt00jIiKM5ZcuXUpERASffvopdevW5cqVK1y5coWEhARbnYIQQogs+u8ABev0dWmXPJHz+qI43b0M81vA4SW2DfBfNk2U4eHhzJgxg4kTJ1K1alV27NjBunXrCA4OBiA2NtakT+UXX3xBSkoKgwYNIiAgwLgMGTLEVqcghBAiix4doOAsJTjU8kco0xJSH8BPA0laMZB9J/+x6Ug+Nu1HaQvSjzL/kr5swlbk2ns6sQn3OB+fRIiPq2HAdL0ednyM2haJBsUxfTCDdUPo36EZ4bVKWO245uYDmw9hJ4QQomAL8HChXinvh7OK2NkRW+0teiSP4roqzLN2F/jJcQzbV83nz4s3czw+SZRCCCFynej4u+zQVyLsQST79WVw19xjtuMMjnzZlxX7TuVoLJIohRBC5DppDX2u4kWX5P8xN+UlALprN1FxXQeunfszx2KRRCmEECLXSWvoYwekYM+UlC50Tx7FNeVOObuLeC1pDn+tyJFYJFEKIYTIlcJrlWDVoPqkDda2U1+ZsAdT2KWviJ0+GQoH5EgckiiFEELkWlWCijDlP11IbmiKcKzpAo42X0pskeo5EoPNx3oVQgghMhNeqwTPl/HlfHwSRy7d4qP1J9ArsFuzhcgOlazaZSQjUqMUQgiR6wV4uBDi42pMkgB6Be+tPJrtgxFIohRCCJEnRMffNSbJNKlKcT7evOmyskoSpRBCiDzhv2PDptFqNIT4ZO/cwpIohRBC5AmPjg2r1WiY3KHiwxF9sok05hFCCJFn/Ldhj3Fs2GwmiVIIIUSeEuDhkiMJMo3cehVCCCEyIYlSCCGEyIQkSiGEECITkiiFEEKITBS4xjxKGXqrJiYm2jgSYW06nY6kpCQSExNllnmRo+Tay5vS8kBaXnicApcob9++DUBQUJCNIxFCCJEb3L59Gw8Pj8e+r1FPSqX5jF6v5/LlyxQuXBiNRpNhmVq1arF///4ci8nax3ua/WVlW0u2Mafsk8o87v3ExESCgoK4ePEi7u7uZsWT2xTkay8r28u1Zz0F8dpTSnH79m2KFSuGnd3jn0QWuBqlnZ0dxYsXz7SMVqvN0Yvd2sd7mv1lZVtLtjGn7JPKPOl9d3f3PPvLqiBfe1nZXq496ymo115mNck00pgnA4MGDcrTx3ua/WVlW0u2Mafsk8rk9PeTkwrytZeV7eXas56Cfu1lpsDdehX5V2JiIh4eHiQkJOTZv+pF3iTXXv4mNUqRbzg5OTFu3DicnJxsHYooYOTay9+kRimEEEJkQmqUQgghRCYkUQohhBCZkEQphBBCZEISpRBCCJEJSZRCCCFEJiRRigIrKSmJ4OBghg8fbutQRAFy+/ZtatWqRdWqValUqRJfffWVrUMST1DghrATIs2HH35InTp1bB2GKGBcXV3Zvn07rq6uJCUlUbFiRTp06IC3t7etQxOPITVKUSCdPn2aEydOEBYWZutQRAGj1WpxdXUF4P79+6Smpj5xmidhW5IoRa6zY8cOXnrpJYoVK4ZGo2H16tXpysyePZvQ0FCcnZ2pUaMGO3futOgYw4cPJzIy0koRi/wkJ66/W7duUaVKFYoXL87IkSPx8fGxUvQiO0iiFLnO3bt3qVKlCp999lmG7y9btoyhQ4cyZswYDh8+TMOGDWnVqhUxMTHGMjVq1KBixYrplsuXL/PTTz9RpkwZypQpk1OnJPKQ7L7+ADw9Pfnzzz+Jjo7mu+++4+rVqzlybiKLlBC5GKBWrVplsq527dqqf//+JuvKlSunRo0aZdY+R40apYoXL66Cg4OVt7e3cnd3VxMmTLBWyCIfyY7r71H9+/dXy5cvz2qIIgdIjVLkKcnJyRw8eJDmzZubrG/evDl79uwxax+RkZFcvHiR8+fP88knn9CvXz/Gjh2bHeGKfMYa19/Vq1dJTEwEDLOO7Nixg7Jly1o9VmE90upV5Cnx8fGkpqbi7+9vst7f358rV67YKCpRUFjj+vvnn3/o06cPSimUUgwePJjKlStnR7jCSiRRijxJo9GYvFZKpVtnjp49e1opIlGQPM31V6NGDf74449siEpkF7n1KvIUHx8ftFptur/e4+Li0v2VL4S1yfVXMEmiFHmKo6MjNWrUICoqymR9VFQU9evXt1FUoqCQ669gkluvIte5c+cOZ86cMb6Ojo7mjz/+wMvLixIlSjBs2DC6d+9OzZo1qVevHl9++SUxMTH079/fhlGL/EKuP5GObRvdCpHe1q1bFZBu6dGjh7HM559/roKDg5Wjo6OqXr262r59u+0CFvmKXH/iURqlZOwkIYQQ4nHkGaUQQgiRCUmUQgghRCYkUQohhBCZkEQphBBCZEISpRBCCJEJSZRCCCFEJiRRCiGEEJmQRCmEEEJkQhKlELnctm3b0Gg03Lp1K8ePrdFo0Gg0eHp6Zlpu/PjxVK1a1eR12rYzZszI1hiFyG6SKIXIRRo3bszQoUNN1tWvX5/Y2Fg8PDxsEtM333zDqVOnLNpm+PDhxMbGUrx48WyKSoicI4OiC5HLOTo6UrRoUZsd39PTEz8/P4u2cXNzw83NDa1Wm01RCZFzpEYpRC7Rs2dPtm/fzsyZM423Lc+fP5/u1uuCBQvw9PRk7dq1lC1bFldXVzp27Mjdu3dZuHAhISEhFClShDfffJPU1FTj/pOTkxk5ciSBgYEUKlSIOnXqsG3btizFOmXKFPz9/SlcuDB9+vTh/v37VvgEhMidJFEKkUvMnDmTevXq0a9fP2JjY4mNjSUoKCjDsklJSfzf//0f33//PRs2bGDbtm106NCBdevWsW7dOhYvXsyXX37JihUrjNv06tWL3bt38/3333PkyBE6depEy5YtOX36tEVxLl++nHHjxvHhhx9y4MABAgICmD179lOduxC5mdx6FSKX8PDwwNHREVdX1yfeatXpdMyZM4dSpUoB0LFjRxYvXszVq1dxc3OjQoUKNGnShK1btxIeHs7Zs2dZunQp//zzD8WKFQMMzxE3bNjAN998w+TJk82Oc8aMGfTu3Zu+ffsCMGnSJDZt2iS1SpFvSY1SiDzI1dXVmCQB/P39CQkJwc3NzWRdXFwcAIcOHUIpRZkyZYzPD93c3Ni+fTtnz5616NjHjx+nXr16JusefS1EfiI1SiHyIAcHB5PXGo0mw3V6vR4AvV6PVqvl4MGD6RrY/De5CiHSk0QpRC7i6Oho0gDHWqpVq0ZqaipxcXE0bNjwqfZVvnx59u3bR0REhHHdvn37njZEIXItSZRC5CIhISH89ttvnD9/Hjc3N7y8vKyy3zJlyvDaa68RERHBp59+SrVq1YiPj2fLli1UqlSJsLAws/c1ZMgQevToQc2aNXnuuedYsmQJx44do2TJklaJVYjcRp5RCpGLDB8+HK1WS4UKFfD19SUmJsZq+/7mm2+IiIjgnXfeoWzZsrRt25bffvvtsS1rHyc8PJyxY8fy7rvvUqNGDS5cuMCAAQOsFqcQuY1GKaVsHYQQInfSaDSsWrWKdu3aZWn7kJAQhg4dmm60ISHyEqlRCiEy1aVLF4uHops8eTJubm5WrRELYStSoxRCPNaZM2cA0Gq1hIaGmr3djRs3uHHjBgC+vr42G6dWCGuQRCmEEEJkQm69CiGEEJmQRCmEEEJkQhKlEEIIkQlJlEIIIUQmJFEKIYQQmZBEKYQQQmRCEqUQQgiRCUmUQgghRCYkUQohhBCZ+H/xB4WRejiIDAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "tm = np.logspace(np.log10(to[0]), np.log10(to[-1]), 100)\n", "hm_0 = w_0.headinside(tm)\n", "plt.semilogx(to, ho / H0, \".\", label=\"obs\")\n", "plt.semilogx(tm, hm_0[0] / H0, label=\"ttim\")\n", "plt.xlabel(\"time [d]\")\n", "plt.ylabel(\"Normalized head (h/H0)\")\n", "plt.title(\"Model results - three layers model\")\n", "plt.legend()\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 7. Create Second Model - multi-layer model\n", "\n", "To investigate whether we can improve the model performance, we will create a multi-layer model. For this, we divide the previous second and third layers into 0.5 m thick layers:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Determine elevation of each layer.\n", "# Thickness of each layer is set to be 0.5 m.\n", "z0 = np.arange(zt, zb, -0.5)\n", "z1 = np.arange(zb, b, -0.5)\n", "zlay = np.append(z0, z1)\n", "zlay = np.append(zlay, b)\n", "zlay = np.insert(zlay, 0, 0)\n", "nlay = len(zlay) - 1 # number of layers\n", "Saq_1 = 1e-4 * np.ones(nlay)\n", "Saq_1[0] = 0.1" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "self.neq 8\n", "solution complete\n" ] } ], "source": [ "ml_1 = ttm.Model3D(\n", " kaq=10, z=zlay, Saq=Saq_1, kzoverkh=1, tmin=1e-5, tmax=0.01, phreatictop=True\n", ")\n", "w_1 = ttm.Well(\n", " ml_1,\n", " xw=0,\n", " yw=0,\n", " rw=rw,\n", " tsandQ=[(0, -Q)],\n", " layers=[1, 2, 3, 4, 5, 6, 7, 8],\n", " rc=rc,\n", " wbstype=\"slug\",\n", ")\n", "ml_1.solve()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 8. Calibration of multi-layer model" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "...................................\n", "Fit succeeded.\n", "[[Fit Statistics]]\n", " # fitting method = leastsq\n", " # function evals = 32\n", " # data points = 27\n", " # variables = 2\n", " chi-square = 0.00108524\n", " reduced chi-square = 4.3410e-05\n", " Akaike info crit = -269.288333\n", " Bayesian info crit = -266.696659\n", "[[Variables]]\n", " kaq0_21: 0.49534703 +/- 0.02256920 (4.56%) (init = 10)\n", " Saq0_21: 4.0607e-04 +/- 1.0398e-04 (25.61%) (init = 0.0001)\n", "[[Correlations]] (unreported correlations are < 0.100)\n", " C(kaq0_21, Saq0_21) = -0.9593\n" ] } ], "source": [ "ca_1 = ttm.Calibrate(ml_1)\n", "ca_1.set_parameter(name=\"kaq0_21\", initial=10, pmin=0)\n", "ca_1.set_parameter(name=\"Saq0_21\", initial=1e-4, pmin=0)\n", "ca_1.seriesinwell(name=\"obs\", element=w_1, t=to, h=ho)\n", "ca_1.fit(report=True)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
optimalstdperc_stdpminpmaxinitialparray
kaq0_210.4953470.0225694.5562410inf10.0000[0.4953470298678706, 0.4953470298678706, 0.495...
Saq0_210.0004060.00010425.6076190inf0.0001[0.00040606934369180614, 0.0004060693436918061...
\n", "
" ], "text/plain": [ " optimal std perc_std pmin pmax initial \\\n", "kaq0_21 0.495347 0.022569 4.556241 0 inf 10.0000 \n", "Saq0_21 0.000406 0.000104 25.607619 0 inf 0.0001 \n", "\n", " parray \n", "kaq0_21 [0.4953470298678706, 0.4953470298678706, 0.495... \n", "Saq0_21 [0.00040606934369180614, 0.0004060693436918061... " ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "RMSE: 0.006339882457078611\n" ] } ], "source": [ "display(ca_1.parameters)\n", "print(\"RMSE:\", ca_1.rmse())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "RMSE has just slightly improved, and the parameter values are more or less similar to the previous values. However, AIC has improved significantly." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 10. Analysis and comparison of simulated values\n", "\n", "We now compare the values in TTim and add the results of the AQTESOLV modelling reported by Yang (2020)." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Comparison of parameter values and error under different models
 k [m/d]Ss [1/m]RMSE
AQTESOLV2.6160000.0000790.001197
ttim-three0.5998110.0002090.006514
ttim-multi0.4953470.0004060.006340
\n" ], "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ta = pd.DataFrame(\n", " columns=[\"k [m/d]\", \"Ss [1/m]\"],\n", " index=[\"AQTESOLV\", \"ttim-three\", \"ttim-multi\"],\n", ")\n", "ta.loc[\"ttim-three\"] = ca_0.parameters[\"optimal\"].values\n", "ta.loc[\"ttim-multi\"] = ca_1.parameters[\"optimal\"].values\n", "ta.loc[\"AQTESOLV\"] = [2.616, 7.894e-5]\n", "ta[\"RMSE\"] = [\n", " 0.001197,\n", " round(ca_0.rmse(), 6),\n", " round(ca_1.rmse(), 6),\n", "]\n", "ta.style.set_caption(\"Comparison of parameter values and error under different models\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "AQTESOLV parameters are quite different from the set parameters in TTim. It also has a better RMSE performance. All TTim models are very similar to each other. However, the multi-layer models performed better." ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## References\n", "\n", "* Batu, V., 1998. Aquifer hydraulics: a comprehensive guide to hydrogeologic data analysis. John Wiley & Sons\n", "* Hyder, Z., Butler Jr, J.J., McElwee, C.D., Liu, W., 1994. Slug tests in partially penetrating wells. Water Resources Research 30, 2945–2957.\n", "* Duffield, G.M., 2007. AQTESOLV for Windows Version 4.5 User's Guide, HydroSOLVE, Inc., Reston, VA.\n", "* Yang, Xinzhu (2020) Application and comparison of different methodsfor aquifer test analysis using TTim. Master Thesis, Delft University of Technology (TUDelft), Delft, The Netherlands." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.4" } }, "nbformat": 4, "nbformat_minor": 4 }