{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 1. Slug Test - Pratt County\n", "**This test is taken from AQTESOLV examples.**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 1. Load 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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this notebook, we reproduce the work of Yang (2020) to check the TTim performance in analysing slug tests. 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 conducted in Pratt County Monitoring Site, US, and reported by Butler (1998). A partially penetrating well is screened in unconsolidated alluvial deposits, consisting of sand and gravel interbedded by clay. The total thickness of the aquifer is 47.87 m. The screen is located at 16.77 m depth and has a screen length 1.52 m the well radius is 0.125 m, and the casing radius 0.064 m.\n", "\n", "The slug displacement is 0.671 m. Head change has been recorded at the slug well.\n", "\n", "The conceptual model can be seen in the figure below." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "jupyter": { "source_hidden": true }, "tags": [ "hide-input" ] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdsAAAE6CAYAAABeYvgTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJSUlEQVR4nO3dd1gU1/s28Htpu7AUQbpS7YoVjKKJir0RNV9bsIA1SOyxhBQFo2LXaGLs5ZdYYyyJLZCoqBENVoxGExVEBSyooKggcN4/fJm47oK7uCuC9+e69rqYM2fOPHPY3Wdn5syMTAghQERERAZjVNIBEBERlXVMtkRERAbGZEtERGRgTLZEREQGxmRLRERkYEy2REREBsZkS0REZGBMtkRERAbGZEtERGRgTLavUUJCAgYMGAAvLy8oFApYWlqiQYMGmDVrFu7evVvS4eldSkoKIiIicPr06deyvhYtWqBFixZa1ZPJZPD29oamG6gdPHgQMpkMMpkMa9as0Vt8a9asgUwmQ1JSks7LRkREQCaT6S0WAEhKSpK2UyaTwcjICOXLl0fHjh0RFxen13UV9V7YvXs3IiIidG7zl19+QWBgIJycnGBmZgY7Ozu0atUK69atw9OnT189aD2YPn06tm/frtc2X/y/vfgqTl++iUJCQuDp6VnSYegNk+1rsnz5cvj6+iI+Ph7jx4/H3r17sW3bNvTo0QNLlizBoEGDSjpEvUtJSUFkZORrS7a6sLKyQmJiIvbt26c2b9WqVbC2ti6BqErGiBEjEBcXh0OHDiEqKgpnzpxBQEAATp06pbd1FPVe2L17NyIjI7VuSwiBAQMG4P3330d+fj7mzZuH3377DWvXrkXdunURFhaGxYsX6y32V2GIZFug4P/24mvw4MEGWR+9GpOSDuBtEBcXh2HDhqFNmzbYvn075HK5NK9Nmzb45JNPsHfv3hKM8O3j7u4OKysrrFq1Cq1atZLKHzx4gB9//BF9+vTB8uXLSzDC18fd3R2NGzcGADRt2hSVK1dGq1atsHjx4kL74PHjx1AoFHrf29bG7NmzsWbNGkRGRmLSpEkq8wIDAzFhwgRcunTptcf1uj3/f6M3H/dsX4Pp06dDJpNh2bJlKom2gJmZGd5//31pOj8/H7NmzUL16tUhl8vh6OiI/v374/r16yrLtWjRAj4+PoiPj8d7770HCwsLeHt7Y8aMGcjPz1epe//+fXzyySfw9vaW2uzYsSMuXLgg1cnJycHUqVOl9To4OGDAgAG4ffu2Sluenp7o3Lkztm3bhjp16kChUMDb2xsLFy6U6hw4cAANGzYEAAwYMEDtEFdhh3w1HTqKjIxEo0aNYGdnB2trazRo0AArV67UeAhYFwMHDsTWrVtx//59qWzjxo0AgN69e2tc5vDhw2jVqhWsrKxgYWGBJk2aYNeuXWr1jh49iqZNm0KhUMDV1RXh4eGFHtrctGkT/P39oVQqYWlpiXbt2ul1r1JXBV/gV69eBfDf4e/o6GgMHDgQDg4OsLCwQHZ2Ni5duoQBAwagSpUqsLCwQIUKFRAYGIizZ89K7RX1XggJCcG3334LACqHQgs71P706VPMnDkT1atXx5dffqmxjrOzM959911p+u7duwgLC0OFChVgZmYGb29vfP7558jOzpbqFBya1XTa4MVDswWH9M+dO4cPP/wQNjY2cHJywsCBA5GRkaGyXFZWFtauXSttV4sWLZCUlAQTExNERUWpravgFMaPP/6ocdt08e+//8La2ho9evRQKd+3bx+MjY1V+m/Tpk1o27YtXFxcYG5ujho1auDTTz9FVlaWyrIhISGwtLTEhQsX0K5dOyiVSri4uGDGjBkAnr3v3333XSiVSlStWhVr165VWb7gvRQTE4MBAwbAzs4OSqUSgYGBuHLlyku3SQiBxYsXo169ejA3N4etrS26d++u1bIlTpBB5ebmCgsLC9GoUSOtlxk6dKgAIIYPHy727t0rlixZIhwcHISbm5u4ffu2VK958+aifPnyokqVKmLJkiUiJiZGhIWFCQBi7dq1Ur3MzExRq1YtoVQqxZQpU8Svv/4qfvrpJzFq1Cixb98+IYQQeXl5on379kKpVIrIyEgRExMjVqxYISpUqCBq1qwpHj16JLXn4eEhKlSoINzd3cWqVavE7t27RZ8+fQQAMXv2bCGEEBkZGWL16tUCgPjiiy9EXFyciIuLE9euXZNib968udq2BwcHCw8PD5WykJAQsXLlShETEyNiYmLEV199JczNzUVkZKRKvcLafFHz5s1FrVq1RGZmplAqlWLx4sXSvEaNGon+/fuL+Ph4AUCsXr1amnfgwAFhamoqfH19xaZNm8T27dtF27ZthUwmExs3bpTqnTt3TlhYWIiaNWuKDRs2iB07doh27doJd3d3AUAkJiZKdadNmyZkMpkYOHCg2Llzp9i6davw9/cXSqVSnDt3Tqo3efJkoe+Pa2Jiosr/rMCZM2cEABEUFCSEENL/sUKFCmLo0KFiz549YsuWLSI3N1fExsaKTz75RGzZskXExsaKbdu2ia5duwpzc3Nx4cIFIUTR74VLly6J7t27CwBSeVxcnHjy5InGmI8cOSIAiIkTJ2q1jY8fPxZ16tQRSqVSzJkzR0RHR4svv/xSmJiYiI4dO6r1xfP/7wIAxOTJk6Xpgv9FtWrVxKRJk0RMTIyYN2+ekMvlYsCAAVK9uLg4YW5uLjp27ChtV8H/tFu3bsLd3V3k5uaqrKtHjx7C1dVVPH36tNBtKoh15syZ4unTp2qv523cuFEAEF9//bUQQojU1FTh5OQkmjdvrrLur776SsyfP1/s2rVLHDhwQCxZskR4eXmJgIAAlfaCg4OFmZmZqFGjhvj6669FTEyMGDBggAAgwsPDRdWqVcXKlSvFr7/+Kjp37iwAiOPHj0vLF7wP3NzcxMCBA8WePXvEsmXLhKOjo3BzcxP37t1TWdeL3wVDhgwRpqam4pNPPhF79+4V69evF9WrVxdOTk4iLS2t0D57EzDZGlhaWpoAIHr37q1V/b///lsAEGFhYSrlx44dEwDEZ599JpU1b95cABDHjh1TqVuzZk3Rrl07aXrKlCkCgIiJiSl0vRs2bBAAxE8//aRSXpB0nk9IHh4eQiaTidOnT6vUbdOmjbC2thZZWVkqy2r6AtMl2T4vLy9PPH36VEyZMkWUL19e5Ofnv7RNTeuuVauWtD4/Pz8hxLMkCUAcOHBAY+yNGzcWjo6O4sGDB1JZbm6u8PHxERUrVpRi6dWrlzA3N1f58Ofm5orq1aurJNvk5GRhYmIiRowYoRLfgwcPhLOzs+jZs6dUZshkW/Cl/eTJE3HixAnRsGFDAUDs2rVLCPHfF2T//v1f2mZubq7IyckRVapUEWPGjJHKi3ovfPzxx1pvW0HyWLJkiVb1lyxZIgCIzZs3q5TPnDlTABDR0dFCiOIl21mzZqnUCwsLEwqFQuU9qVQqRXBwsFqb+/fvFwDEtm3bpLIbN24IExMTtR+RLyqItbDXoUOHVOoPGzZMmJmZibi4ONGyZUvh6OgoUlJSCm0/Pz9fPH36VMTGxgoA4syZM9K84OBgte+Jp0+fCgcHBwFAnDx5UipPT08XxsbGYuzYsVJZwXupW7duKuv8448/BAAxdepUlXU9/10QFxcnAIi5c+eqLHvt2jVhbm4uJkyYUGS/lTQeRn7D7N+/H8CzwzXPe+edd1CjRg38/vvvKuXOzs545513VMrq1KkjHQIEgD179qBq1apo3bp1oevduXMnypUrh8DAQOTm5kqvevXqwdnZGQcOHFCpX6tWLdStW1elLCgoCJmZmTh58qS2m6uVffv2oXXr1rCxsYGxsTFMTU0xadIkpKen49atW6/U9sCBA3H8+HGcPXsWK1euRKVKldCsWTO1ellZWTh27Bi6d+8OS0tLqdzY2Bj9+vXD9evXcfHiRQDP/oetWrWCk5OTSr1evXqptPnrr78iNzcX/fv3V+lzhUKB5s2bq/X5ywghVNrJzc3VarmJEyfC1NQUCoUCvr6+SE5OxtKlS9GxY0eVev/73//Uls3NzcX06dNRs2ZNmJmZwcTEBGZmZvj333/x999/6xS/Iezbtw9KpRLdu3dXKS/4fL34edLF86d+gGefuydPnmj1nmzRogXq1q0rHUIHgCVLlkAmk2Ho0KFarX/UqFGIj49Xe9WrV0+l3vz581GrVi0EBATgwIED+OGHH+Di4qJS58qVKwgKCoKzs7P0GWvevDkAqP0fZTKZynvDxMQElStXhouLC+rXry+V29nZwdHRUeW7qECfPn1Upps0aQIPDw/p+0+TnTt3QiaToW/fvirvcWdnZ9StW1fnz8vrxgFSBmZvbw8LCwskJiZqVT89PR0A1D4MAODq6qr2xi1fvrxaPblcjsePH0vTt2/fhru7e5HrvXnzJu7fvw8zMzON8+/cuaMy7ezsrFanoKxgG/Thzz//RNu2bdGiRQssX74cFStWhJmZGbZv345p06apbGdxNGvWDFWqVMHSpUuxefNmjB49WuOgn3v37kEIUej/Bfhvu9PT04vsnwI3b94EAOl85ouMjHT7LRwbG4uAgACVssTExJdePjFq1Cj07dsXRkZGKFeuHLy8vDT2gaZtHzt2LL799ltMnDgRzZs3h62tLYyMjDB48OBX/t9oUvA+1uXz5OzsrLY9jo6OMDExeaX36oufvYLxGNpu98iRIzF48GBcvHgR3t7eWL58Obp3767xvaNJxYoV4efn99J6crkcQUFBGD9+PBo0aIA2bdqozH/48CHee+89KBQKTJ06FVWrVoWFhQWuXbuGDz74QG17LCwsoFAoVMoKLr16kZmZGZ48eaJWXtjno6j/x82bNyGEUPkR+zxvb+9Cl30TMNkamLGxMVq1aoU9e/bg+vXrqFixYpH1Cz7AqampanVTUlJgb2+vcwwODg5qg6teZG9vj/Llyxc6KtrKykplOi0tTa1OQZmmHwAvUigUKoNJCryY1Ddu3AhTU1Ps3LlT5QOuz8spBgwYgC+++AIymQzBwcEa6xQkkdTUVLV5KSkpACD9b8qXL19k/xQoqL9lyxZ4eHi80jYAkC4te17BD4GiaPulrSkB//DDD+jfvz+mT5+uUn7nzh2UK1fupW3qys/PD3Z2dtixYweioqJeOhq6fPnyOHbsGIQQKnVv3bqF3Nxc6X9Q8N56ftAUoN8fji8KCgrCxIkT8e2336Jx48ZIS0vDxx9/rPf1/PXXX5g0aRIaNmyI+Ph4zJs3D2PHjpXm79u3DykpKThw4IC0NwtAZeCgvhX2+ahcuXKhy9jb20Mmk+HQoUMaB5pqKnuT8DDyaxAeHg4hBIYMGYKcnBy1+U+fPsUvv/wCAGjZsiWAZ19iz4uPj8fff/+tcpmKtjp06IB//vlH4zWlBTp37oz09HTk5eXBz89P7VWtWjWV+ufOncOZM2dUytavXw8rKys0aNAAQNG/9D09PfHPP/+ofLmlp6fjyJEjKvVkMhlMTExgbGwslT1+/Bjff/+9llv/csHBwQgMDMT48eNRoUIFjXWUSiUaNWqErVu3qmxPfn4+fvjhB1SsWBFVq1YFAAQEBOD333+X9lwBIC8vD5s2bVJps127djAxMcHly5c19rk2CfB5VlZWassXdqRCX2QymdqX3K5du3Djxg2VsqLeC7rsEZqammLixIm4cOECvvrqK411bt26hT/++AMA0KpVKzx8+FDtx9n//d//SfMBwMnJCQqFAgkJCSr1duzY8dKYivLiUabnKRQKDB06FGvXrsW8efNQr149NG3a9JXW96KsrCz06NEDnp6e2L9/P4YPH45PP/0Ux44dk+oU/Ah58f+4dOlSvcbyvHXr1qlMHzlyBFevXi3ypjSdO3eGEAI3btzQ+FmpXbu2weLVB+7Zvgb+/v747rvvEBYWBl9fXwwbNgy1atXC06dPcerUKSxbtgw+Pj4IDAxEtWrVMHToUCxatAhGRkbo0KEDkpKS8OWXX8LNzQ1jxozRef2jR4/Gpk2b0KVLF3z66ad455138PjxY8TGxqJz584ICAhA7969sW7dOnTs2BGjRo3CO++8A1NTU1y/fh379+9Hly5d0K1bN6lNV1dXvP/++4iIiICLiwt++OEHxMTEYObMmbCwsAAAVKpUCebm5li3bh1q1KgBS0tLuLq6wtXVFf369cPSpUvRt29fDBkyBOnp6Zg1a5bazSQ6deqEefPmISgoCEOHDkV6ejrmzJmj11+xrq6uWu0pR0VFoU2bNggICMC4ceNgZmaGxYsX46+//sKGDRukL60vvvgCP//8M1q2bIlJkybBwsIC3377rdplFJ6enpgyZQo+//xzXLlyBe3bt4etrS1u3ryJP//8E0qlUqebPZSEzp07Y82aNahevTrq1KmDEydOYPbs2WpHZYp6LxR8Sc6cORMdOnSAsbEx6tSpU+gPhfHjx+Pvv//G5MmT8eeffyIoKAhubm7IyMjAwYMHsWzZMkRGRqJp06bo378/vv32WwQHByMpKQm1a9fG4cOHMX36dHTs2FEax1BwLnDVqlWoVKkS6tatiz///BPr169/pf6pXbs2Dhw4gF9++QUuLi6wsrJS+eEaFhaGWbNm4cSJE1ixYoVObScnJ+Po0aNq5Q4ODqhUqRIAIDQ0FMnJydL7ae7cuYiLi0Pv3r1x6tQplCtXDk2aNIGtrS1CQ0MxefJkmJqaYt26dWo/pvXp+PHjGDx4MHr06IFr167h888/R4UKFRAWFlboMk2bNsXQoUMxYMAAHD9+HM2aNYNSqURqaioOHz6M2rVrY9iwYQaL+ZWV5Oist83p06dFcHCwcHd3F2ZmZkKpVIr69euLSZMmiVu3bkn18vLyxMyZM0XVqlWFqampsLe3F3379pUumynw/Kja52ka0Xvv3j0xatQo4e7uLkxNTYWjo6Po1KmTdHmGEM9GFc6ZM0fUrVtXKBQKYWlpKapXry4++ugj8e+//0r1PDw8RKdOncSWLVtErVq1hJmZmfD09BTz5s1Ti2XDhg2ievXqwtTUVG1U59q1a0WNGjWEQqEQNWvWFJs2bdIY+6pVq0S1atWEXC4X3t7eIioqSqxcuVLtMprijEYuTGGjZw8dOiRatmwplEqlMDc3F40bNxa//PKL2vJ//PGHaNy4sZDL5cLZ2VmMHz9eLFu2TC1mIYTYvn27CAgIENbW1kIulwsPDw/RvXt38dtvv0l1XuelPy8qGEEaHx+vNu/evXti0KBBwtHRUVhYWIh3331XHDp0SOP/orD3QnZ2thg8eLBwcHAQMplMYx9psmPHDtGpUyfh4OAgTExMhK2trQgICBBLliwR2dnZUr309HQRGhoqXFxchImJifDw8BDh4eFqlxdlZGSIwYMHCycnJ6FUKkVgYKBISkoqdDTy85fhPd9Pz8d++vRp0bRpU2FhYSEAaHx/tmjRQtjZ2alcXleUl41G7tOnjxBCiOXLl2t8D1+6dElYW1uLrl27SmVHjhwR/v7+wsLCQjg4OIjBgweLkydPqi0fHBwslEqlWkyFfaYKvite7KPo6GjRr18/Ua5cOenyqOe/YwrWpenKhFWrVolGjRpJn8FKlSqJ/v37q1xi9CaSCfGKdwagt46npyd8fHywc+fOkg6FqFS7desWPDw8MGLECMyaNaukwzG4NWvWYMCAAYiPj9f5NElpx8PIRESv2fXr13HlyhXMnj0bRkZGGDVqVEmHRAbGAVJERK/ZihUr0KJFC5w7dw7r1q0rdGAelR08jExERGRgpWbPNioqCg0bNoSVlRUcHR3RtWtX6Y49BYQQiIiIgKurK8zNzaVfjkRERCWp1CTb2NhYfPzxxzh69ChiYmKQm5uLtm3bqlxOMWvWLMybNw/ffPMN4uPj4ezsjDZt2uDBgwclGDkREb3tSu1h5Nu3b8PR0RGxsbFo1qwZhBBwdXXF6NGjMXHiRADP7gbj5OSEmTNn4qOPPirhiImI6G1VakcjF9zqr+B+nImJiUhLS0Pbtm2lOnK5HM2bN8eRI0cKTbbZ2dkqdzHKz8/H3bt3Ub58+RJ5MDYREZU8IQQePHgAV1dXne9TrkmpTLZCCIwdOxbvvvsufHx8APx3r80Xb1Lt5OSk8akTBaKiot74u/QQEVHJuHbt2kvvaa+NUplshw8fjoSEBBw+fFht3ot7o+KFG5C/KDw8XOWm3BkZGXB3d8fx49dgaWld6HJERFR2PXyYCT8/N7WHsBRXqUu2I0aMwM8//4yDBw+q/NooeGRTWlqayqPAbt26VegjmYBnh5o13WfX0tIaVlZMtkREbzN9nU4sNaORhRAYPnw4tm7din379sHLy0tlvpeXF5ydnRETEyOV5eTkIDY2Fk2aNHnd4RIREUlKzZ7txx9/jPXr12PHjh2wsrKSztHa2NjA3NwcMpkMo0ePxvTp01GlShVUqVIF06dPh4WFBYKCgko4eiIiepuVmmT73XffAYDa8w5Xr16NkJAQAMCECRPw+PFjhIWF4d69e2jUqBGio6P1dsydiIioOErtdbaGkpmZCRsbG1y4kMFztkREb6kHDzJRvboNMjIy1J6zXRyl5pwtERFRacVkS0REZGBMtkRERAbGZEtERGRgTLZEREQGxmRLRERkYEy2REREBsZkS0REZGBMtkRERAbGZEtERGRgTLZEREQGxmRLRERkYEy2RG+o33/fjQoVZNLL3d0EtWqVxwcfNMOKFV8jOzvb4DFkZT3EpEmj0aCBK7y9FWjTph527NioUxt//nkY/fp1RM2atqhUyRxNm1bB/PlfSfNHjw5R2c4XXydOHJXqPnz4AFOnTsCHH7ZF7doOqFBBhrlzI/S1uUQGU2oesUf0tjl79iQAYMWKrXB0dEFeXi5u376J2NhofPXVOGzevAY//rgfNjblDBbD4MEf4MyZeISHz4C3d1Vs374eYWEfIj8/H926vfw50du2rcfIkf0QGNgTCxb8H5RKS1y9ehk3b6ZIdUaP/hL9+oWqLRsSEggzMznq1Wsold27l45165ahZs26aN++K9avX6GfDSUyMCZbojfUX3+dhFJpifbtu0Imk0nlnTr9D35+TTBmTAjmzJmEr75aaJD1//77bhw8GINvv12Prl0/BAA0bRqA69evYurU8Xj//V4wNjYudPnU1BuYMGEo+vb9CFFRi6Xypk0DVOp5elaCp2cllbK4uFjcvXsHo0Z9obKOihU9cP78PchkMty9e4fJlkoNHkYmKiH5+flFzk9IOIFateqpJNoCPXsGw8WlInbv/slQ4WHv3m1QKi3RuXMPlfJevQYgLS0FJ08eK3L5DRtW4NGjLHz88USd171hw0rIZDL07j1QpVwmk2nsD23NnRuBChVkOH8+AUOH9kD16jaoVcsOERFjkZubi0uXLqJPn/aoWtUKjRp5YvHiWcVeF9HzmGyJSkBiUhpatQnHzl1/apx/9246btxIho9P/ULb8Pauips3UzUmbSEEcnNztXoV5sKFv1ClSg2YmKgeAKtRow4A4OLFv4rcxqNHD6JcOTtcunQBbdrUg7u7CerUccTEiaF48CCz0OUyMzOwa9cWvPtuK7i7exW5juIKDe2JmjXrYtmynxAUNATLl89HRMQYDBrUFa1adcKKFdvQtGlLTJs2Ebt3bzVIDPR24WFkotfs9JnL6NN3FjIyHyE0bCEiJ/fDoIHtVOr89dez87VFJducnGyYm1vAyEj9N3NcXCx69AjQsJS6o0cT4ebmqVZ+7146PDy81crLlbOT5hclLe0Gnjx5hI8+6oHhw8Ph67sAZ87EY86cybh48S9s23ZI417q9u0b8OTJY3z44SCt4i+OPn2G4qOPxgIAmjVrjdjYaKxe/Q1WrNiKDh26AQCaNGmB337biW3b1qFjxw8MFgu9HZhsiV6jx09y0ClwskrZpIjv4e3tjIAWdaWyhIQTAIBateoV2lZS0iVUqlRN47w6dXyxe3e8VjE5ObkWOq+oQ7YvO5ybn5+PJ0+eIDx8MoYP/xTAswRmamqGyZNH49Ch39GsWWu15TZuXAlb2/Jo376bVvEXR+vWnVWmq1SpgfPnzyAgoINUZmJiAk/Pyrh+/arB4qC3B5Mt0WtkrjBDUO8W2LL1MHJynh3CrVHDHQ3qV1apd/bsSZiamqJq1Voa2zl9Oh63b99Enz5DNc5XKi2LTNTPe/EwcQFb2/Ia917v378L4L893MLY2pZHYuK/aNFCda89IKADJk8ejb/+OqmWbM+fT8CZM8cxaNAoyOVyreIvDltb1dhNTc1gbm4BhUKhUm5mZoaHDws/5E2kLZ6zJXrNZs8aDHd3R2l6zcqxsLFRqtT566+TqFq1FszMzNSWF0Jg1qwvoFCYIzg4TOM64uJi4eFhqtXr2rUkjW3UqFEb//77t9p53QsXzgIAqlXzKXI7C87taoofgMbD3xs3rgQABAUNLrJtotKGe7ZEb5jMzAxcvXoFPXuGqM17+vQpvvhiBGJjozFjxhI4OjprbEMfh5Hbt++GdeuWY9eun9ClSy+p/Mcf18LZ2RUNGjQqst1Onf6HdeuWYd++PSrnnvft2w0AaNCgsUr97Oxs/PTTD6hf/x1Ur150IicqbZhsid4wZ8+ehBAClpZWOHHiKPLz85GRcQ8JCcexefMa3LqViunTF6Nfv48KbcPS0gp16/q9UhwtW3ZAs2Zt8Nlnw/DwYSY8PStjx44N2L9/LxYt+kG6/jUuLha9erXCmDGTMGbMJGn55s3bok2bQCxYMAX5+flo0KAxEhKOY/78SLRu3RnvvPOuyvr27t2O+/fv4sMPi96r3bdvDx49ykJW1gMAwD//nMfOnVsAAK1adYS5ucUrbTeRIZSqZHvw4EHMnj0bJ06cQGpqKrZt24auXbtK84UQiIyMxLJly3Dv3j00atQI3377LWrV0nzei+hNVHDnqJUrF2LlyoWQy+WwsbFFpUrV0bNnCPr0GQonJ5fXEsuKFVsxc+bnmDNnEu7fv4tKlapj8eIN6NKlt1RHCIG8vDyNlyB9990mzJ8fiXXrlmH+/Eg4Obli8OAxGDt2slrdjRtXwsJCqdK2JuHhw1QGLe3c+SN27vwRQOEjq4lKmkwUnEApBfbs2YM//vgDDRo0wP/+9z+1ZDtz5kxMmzYNa9asQdWqVTF16lQcPHgQFy9ehJWVlVbryMzMhI2NDS5cyICVlbWBtoTeds1bTsClS89uWXjsyAJUrGhfwhER0fMePMhE9eo2yMjIgLX1q+eCUrVn26FDB3To0EHjPCEEFixYgM8//xwffPDsmri1a9fCyckJ69evx0cfFX7IjYiIyJBKVbItSmJiItLS0tC2bVupTC6Xo3nz5jhy5EihyTY7O1vl6SmZmRzmT6/uxo1k3L17p9D5T57cBMSzS2guXDiDe/dsC61rZ2ePChXc9R4jEb0+ZSbZpqWlAQCcnJxUyp2cnHD1auEXpUdFRSEyMtKgsdHb5caNZDRvXgOPHz/Sqn5w8J4i55ubWyA29m8mXKJSrMwk2wIv3tVGCFHknW7Cw8MxduxYaTozMxNubm4Gi4/Kvrt37+Dx40eIiIiAp6fnK7WVlJSEiIgI3L17h8mWqBQrM8nW2fnZ9YZpaWlwcflvpOatW7fU9nafJ5fLDXqnGnp7eXp6onr16iUdBhG9AcrMHaS8vLzg7OyMmJgYqSwnJwexsbFo0qRJCUZGRERvu1K1Z/vw4UNcunRJmk5MTMTp06dhZ2cHd3d3jB49GtOnT0eVKlVQpUoVTJ8+HRYWFggKCirBqImI6G1XqpLt8ePHERDw32PDCs61BgcHY82aNZgwYQIeP36MsLAw6aYW0dHRWl9jS0REZAilKtm2aNECRd2DQyaTISIiAhEREa8vKCIiopcoM+dsiYiI3lRMtkRERAbGZEtERGRgTLZEREQGxmRLRERkYEy2REREBsZkS0REZGBMtkRERAbGZEtERGRgTLZEREQGptXtGp9/3qu2vvjiC9jZ2em8HBERUVmjVbJdsGAB/P39YWZmplWjhw8fxvDhw5lsiYiIoMODCLZt2wZHR0et6vIpO0RERP/R6pzt6tWrYWNjo3WjS5cuhZOTU7GDIiIiKku02rMNDg7WqVE+rJ2IiOg/r/Q824cPHyI/P1+lzNra+pUCIiIiKmt0vvQnMTERnTp1glKphI2NDWxtbWFra4ty5crB1tbWEDESERGVajrv2fbp0wcAsGrVKjg5OUEmk+k9KCIiorJE52SbkJCAEydOoFq1aoaIh4iIqMzR+TByw4YNce3aNUPEQkREVCbpvGe7YsUKhIaG4saNG/Dx8YGpqanK/Dp16ugtOCIiorJA52R7+/ZtXL58GQMGDJDKZDIZhBCQyWTIy8vTa4BERESlnc7JduDAgahfvz42bNjAAVJERERa0DnZXr16FT///DMqV65siHiIiIjKHJ0HSLVs2RJnzpwxRCx6s3jxYnh5eUGhUMDX1xeHDh0q6ZCIiOgtpvOebWBgIMaMGYOzZ8+idu3aagOk3n//fb0FVxybNm3C6NGjsXjxYjRt2hRLly5Fhw4dcP78ebi7u5dobERE9HbSOdmGhoYCAKZMmaI2700YIDVv3jwMGjQIgwcPBvDs8YC//vorvvvuO0RFRWndzqNHWTA2NjZUmFSGPXnyGACQnZ2Nx48fv1Jb2dnZUpuPHmW9cmxEpB19f95kQgih1xZLUE5ODiwsLPDjjz+iW7duUvmoUaNw+vRpxMbGqi2TnZ0tfaEBQGZmJtzc3F5LvERE9GbLyMjQyz3/dT5n+ya7c+cO8vLy1B7v5+TkhLS0NI3LREVFwcbGRnox0RIRkb5pdRh54cKFGDp0KBQKhVaNLlmyBH369Cmxh8i/eDlSwTXAmoSHh2Ps2LHSdMGe7f4NX8LSQrvtJXpR6u17yMh8hEfKusVuwyLr2UBEG2sLuDjwIR9Er9PDR08Q8OFXemtPq2Q7ZswYfPjhh1on2wkTJqBt27avPdna29vD2NhYbS/21q1bhT7MXi6XQy6Xq5VbKMxgYW5mkDip7Kvk/uz99sC6+MnWKjNTX+EQkY5efHzsq9Iq2Qoh0KpVK5iYaDee6lUHhRSXmZkZfH19ERMTo3LONiYmBl26dCmRmIiIiLTKnpMnT9ap0S5dusDOzq5YAb2qsWPHol+/fvDz84O/vz+WLVuG5ORkaRQ1ERHR62aQZFuSevXqhfT0dEyZMgWpqanw8fHB7t274eHhUdKhERHRW0rn62xLg7CwMISFhZV0GERERADK2KU/REREbyImWyIiIgNjsiUiIjIwnZPtlClT8OjRI7Xyx48fa7xfMhER0dtO52QbGRmJhw8fqpU/evQIkZGRegmKiIioLNE52RZ268MzZ86U2LW1REREbzKtL/2xtbWFTCaDTCZD1apVVRJuXl4eHj58yBtHEBERaaB1sl2wYAGEEBg4cCAiIyNhY2MjzTMzM4Onpyf8/f0NEiQREVFppnWyDQ4OBgB4eXmhSZMmMDU1NVhQREREZYnOd5Bq3rw58vPz8c8//+DWrVtqT0Zo1qyZ3oIjIiIqC3ROtkePHkVQUBCuXr0KIYTKPJlMhry8PL0FR0REVBbonGxDQ0Ph5+eHXbt2wcXFpdCHshMREdEzOifbf//9F1u2bEHlypUNEQ8REVGZo/N1to0aNcKlS5cMEQsREVGZpNWebUJCgvT3iBEj8MknnyAtLQ21a9dWG5Vcp04d/UZIRERUymmVbOvVqweZTKYyIGrgwIHS3wXzOECKiIhInVbJNjEx0dBxEBERlVlaJVsPDw9Dx0FERFRm6Twa+eeff9ZYLpPJoFAoULlyZXh5eb1yYERERGWFzsm2a9euaudvAdXztu+++y62b98OW1tbvQVKRERUWul86U9MTAwaNmyImJgYZGRkICMjAzExMXjnnXewc+dOHDx4EOnp6Rg3bpwh4iUiIip1dN6zHTVqFJYtW4YmTZpIZa1atYJCocDQoUNx7tw5LFiwQGW0MhER0dtM5z3by5cvw9raWq3c2toaV65cAQBUqVIFd+7cefXoiIiIygCdk62vry/Gjx+P27dvS2W3b9/GhAkT0LBhQwDPbulYsWJF/UUJYNq0aWjSpAksLCxQrlw5jXWSk5MRGBgIpVIJe3t7jBw5Ejk5OXqNg4iISFc6H0ZeuXIlunTpgooVK8LNzQ0ymQzJycnw9vbGjh07AAAPHz7El19+qddAc3Jy0KNHD/j7+2PlypVq8/Py8tCpUyc4ODjg8OHDSE9PR3BwMIQQWLRokV5jISIi0oXOybZatWr4+++/8euvv+Kff/6BEALVq1dHmzZtYGT0bEe5a9eu+o4TkZGRAIA1a9ZonB8dHY3z58/j2rVrcHV1BQDMnTsXISEhmDZtmsZD30RERK+DzskWeHaZT/v27dG+fXt9x1NscXFx8PHxkRItALRr1w7Z2dk4ceIEAgICNC6XnZ2N7OxsaTozM9PgsRIR0dtFq2S7cOFCDB06FAqFAgsXLiyy7siRI/USmK7S0tLg5OSkUmZrawszMzOkpaUVulxUVJS010xERGQIWiXb+fPno0+fPlAoFJg/f36h9WQymU7JNiIi4qWJLj4+Hn5+flq1p+lB9gU32ihMeHg4xo4dK01nZmbCzc1Nq/URERFpQ+cHEejzoQTDhw9H7969i6zj6empVVvOzs44duyYStm9e/fw9OlTtT3e58nlcsjlcq3WQUREVBzFOmcLPBsdnJiYiEqVKsHEpHjN2Nvbw97evrghqPD398e0adOQmpoKFxcXAM8GTcnlcvj6+uplHURERMWh83W2jx49wqBBg2BhYYFatWohOTkZwLNztTNmzNB7gAWSk5Nx+vRpJCcnIy8vD6dPn8bp06fx8OFDAEDbtm1Rs2ZN9OvXD6dOncLvv/+OcePGYciQIRyJTEREJUrnZBseHo4zZ87gwIEDUCgUUnnr1q2xadMmvQb3vEmTJqF+/fqYPHkyHj58iPr166N+/fo4fvw4AMDY2Bi7du2CQqFA06ZN0bNnT3Tt2hVz5swxWExERETa0Pn47/bt27Fp0yY0btxYZeBRzZo1cfnyZb0G97w1a9YUeo1tAXd3d+zcudNgMRARERWHznu2t2/fhqOjo1p5VlZWkaN+iYiI3lY6J9uGDRti165d0nRBgl2+fDn8/f31FxkREVEZofNh5KioKLRv3x7nz59Hbm4uvv76a5w7dw5xcXGIjY01RIxERESlms57tk2aNMEff/yBR48eoVKlSoiOjoaTkxPi4uJ4iQ0REZEGxbpAtnbt2li7dq2+YyEiIiqTipVs8/PzcenSJdy6dQv5+fkq85o1a6aXwIiIiMoKnZPt0aNHERQUhKtXr0IIoTJPJpMhLy9Pb8ERERGVBTon29DQUPj5+WHXrl1wcXHh5T5EREQvoXOy/ffff7FlyxZUrlzZEPEQERGVOTqPRm7UqBEuXbpkiFiIiIjKJK32bBMSEqS/R4wYgU8++QRpaWmoXbs2TE1NVerWqVNHvxESERGVclol23r16kEmk6kMiBo4cKD0d8E8DpAiIiJSp/PD44mIiEg3WiVbDw8PQ8dBRERUZuk8QIqIiIh0U6w7SBHRy/X8+Gvcuj+r0PkbN26EUqmUptevX48NGzZI00YiB/a2Vtj87SiDxklEhsdkS2Qgd+49wO07GVrXz8rKwu3btw0YERGVFCZbIgMzMjJC+fLlX1pPqVTCwcEBAJCenq5233EiKr2KlWzv37+PLVu24PLlyxg/fjzs7Oxw8uRJODk5oUKFCvqOkahUK1++PH755ZeX1gsKCkJQUBAAIDAwkHu5RGWIzsk2ISEBrVu3ho2NDZKSkjBkyBDY2dlh27ZtuHr1Kv7v//7PEHESERGVWjon27FjxyIkJASzZs2ClZWVVN6hQwfpVzkRAYsiQpAh91G7yxoRvX10Trbx8fFYunSpWnmFChWQlpaml6CIyoJaVSti6c6zKiOMC1OtWjXMmTNHmraxsUH2k0eGDI+IXiOdk61CoUBmZqZa+cWLF6XBHUT0TM79K1qde3V2KAez7KvS9IKpo/BZ5HzcvJ1uyPCI6DXROdl26dIFU6ZMwebNmwE8uy9ycnIyPv30U/zvf//Te4BEpZmVAnCyt3lpvfLWppBnJ0vTXnaAnY0CeblWRSxFRKWFTDz/dAEtZGZmomPHjjh37hwePHgAV1dXpKWlwd/fH7t371a5SF9fkpKS8NVXX2Hfvn1IS0uDq6sr+vbti88//xxmZmZSveTkZHz88cfYt28fzM3NERQUhDlz5qjU0Wb7bGxscGzbV7BUKvS+LURE9OZ7mPUEjbp9iYyMDFhbW79yezrv2VpbW+Pw4cPYt28fTp48ifz8fDRo0ACtW7d+5WAKc+HCBeTn52Pp0qWoXLky/vrrLwwZMgRZWVnSea68vDx06tQJDg4OOHz4MNLT0xEcHAwhBBYtWmSw2IiIiF5G5z3bpKQkeHp6Gigc7c2ePRvfffcdrly5AgDYs2cPOnfujGvXrsHV1RXAs9vhhYSE4NatW1r/MuGeLRER6XvPVucHEXh7e+Pdd9/F0qVLcffu3VcOoLgyMjJgZ2cnTcfFxcHHx0dKtADQrl07ZGdn48SJE4W2k52djczMTJUXERGRPumcbI8fPw5/f39MnToVrq6u6NKlC3788UdkZ2cbIj6NLl++jEWLFiE0NFQqS0tLg5OTk0o9W1tbmJmZFXlJUlRUFGxsbKSXm5ubweImIqK3k87JtkGDBpg9ezaSk5OxZ88eODo64qOPPoKjoyMGDhyoU1sRERGQyWRFvo4fP66yTEpKCtq3b48ePXpg8ODBKvNkMpnaOoQQGssLhIeHIyMjQ3pdu3ZNp20gIiJ6mWI/iEAmkyEgIAABAQEYNmwYBg0ahLVr12LVqlVatzF8+HD07t27yDrPnx9OSUlBQEAA/P39sWzZMpV6zs7OOHbsmErZvXv38PTpU7U93ufJ5XLI5XKtYyYiItJVsZPttWvXsGHDBqxfvx5nz56Fv78/vvnmG53asLe3h729vVZ1b9y4gYCAAPj6+mL16tUwMlLdKff398e0adOQmpoKFxcXAEB0dDTkcjl8fX11iouIiEifdE62y5Ytw7p16/DHH3+gWrVq6NOnD7Zv327QEcopKSlo0aIF3N3dMWfOHJU78jg7OwMA2rZti5o1a6Jfv36YPXs27t69i3HjxmHIkCF6GUlGRERUXDon26+++gq9e/fG119/jXr16hkgJHXR0dG4dOkSLl26hIoVK6rMK7hyydjYGLt27UJYWBiaNm2qclMLIiKikqTzdbYvG3BU2vE6WyIiKpE7SCUkJMDHxwdGRkY4e/ZskXXr1KnzykERERGVJVol23r16iEtLQ2Ojo6oV68eZDIZnt8hLpiWyWTIy8szWLBERESlkVbJNjExUXp8XmJiokEDIiIiKmu0SrYeHh7S31evXkWTJk1gYqK6aG5uLo4cOaJSl4iIiIpxB6mAgACN90TOyMhAQECAXoIiIiIqS3ROtoWNRk5PTzfIs2yJiIhKO62vs/3ggw8APBsMFRISonKLw7y8PCQkJKBJkyb6j5CIiKiU0zrZ2tjYAHi2Z2tlZQVzc3NpnpmZGRo3bowhQ4boP0IiIqJSTutku3r1agDPHgwwbtw4HjImIiLSks63a5w8ebIh4iAiIiqzivXUny1btmDz5s1ITk5GTk6OyryTJ0/qJTAiIqKyQufRyAsXLsSAAQPg6OiIU6dO4Z133kH58uVx5coVdOjQwRAxEhERlWo6J9vFixdj2bJl+Oabb2BmZoYJEyYgJiYGI0eOREZGhiFiJCIiKtV0TrbJycnSJT7m5uZ48OABAKBfv37YsGGDfqMjIiIqA3ROts7OzkhPTwfw7DaOR48eBfDsnsk6Pq2PiIjoraBzsm3ZsiV++eUXAMCgQYMwZswYtGnTBr169UK3bt30HiAREVFpp/No5GXLliE/Px8AEBoaCjs7Oxw+fBiBgYEIDQ3Ve4BERESlnc7J1sjICEZG/+0Q9+zZEz179tRrUERERGWJVsk2ISFB6wbr1KlT7GCIiIjKIq2Sbb169SCTyV46AEomkyEvL08vgREREZUVWiXbxMREQ8dBRERUZmmVbD08PAwdBxERUZml86U/APD999+jadOmcHV1xdWrVwEACxYswI4dO/QaHBERUVmgc7L97rvvMHbsWHTs2BH379+XztGWK1cOCxYs0Hd8kvfffx/u7u5QKBRwcXFBv379kJKSolInOTkZgYGBUCqVsLe3x8iRI9UelEBERPS66ZxsFy1ahOXLl+Pzzz+HsbGxVO7n54ezZ8/qNbjnBQQEYPPmzbh48SJ++uknXL58Gd27d5fm5+XloVOnTsjKysLhw4exceNG/PTTT/jkk08MFhMREZE2dL7ONjExEfXr11crl8vlyMrK0ktQmowZM0b628PDA59++im6du2Kp0+fwtTUFNHR0Th//jyuXbsGV1dXAMDcuXMREhKCadOmwdra2mCxERERFUXnPVsvLy+cPn1arXzPnj2oWbOmPmJ6qbt372LdunVo0qQJTE1NAQBxcXHw8fGREi0AtGvXDtnZ2Thx4kShbWVnZyMzM1PlRUREpE86J9vx48fj448/xqZNmyCEwJ9//olp06bhs88+w/jx4w0Ro2TixIlQKpUoX748kpOTVQZkpaWlwcnJSaW+ra0tzMzMkJaWVmibUVFRsLGxkV5ubm4Gi5+IiN5OOifbAQMGYPLkyZgwYQIePXqEoKAgLFmyBF9//TV69+6tU1sRERGQyWRFvo4fPy7VHz9+PE6dOoXo6GgYGxujf//+KjfakMlkausQQmgsLxAeHo6MjAzpde3aNZ22gYiI6GV0PmcLAEOGDMGQIUNw584d5Ofnw9HREQBw48YNVKhQQet2hg8f/tIE7enpKf1tb28Pe3t7VK1aFTVq1ICbmxuOHj0Kf39/ODs749ixYyrL3rt3D0+fPlXb432eXC6HXC7XOmYiIiJdFSvZFrC3twfw7BDutGnTsGLFCjx+/Fin5Qva0FXBHm12djYAwN/fH9OmTUNqaipcXFwAANHR0ZDL5fD19S3WOoiIiPRB68PI9+/fR58+feDg4ABXV1csXLgQ+fn5mDRpEry9vXH06FGsWrXKIEH++eef+Oabb3D69GlcvXoV+/fvR1BQECpVqgR/f38AQNu2bVGzZk3069cPp06dwu+//45x48ZhyJAhHIlMREQlSus9288++wwHDx5EcHAw9u7dizFjxmDv3r148uQJ9uzZg+bNmxssSHNzc2zduhWTJ09GVlYWXFxc0L59e2zcuFE6BGxsbIxdu3YhLCwMTZs2hbm5OYKCgjBnzhyDxUVERKQNmXjZo3z+Pw8PD6xcuRKtW7fGlStXULlyZYwcOdKgd40qCZmZmbCxscGxbV/BUqko6XCIiKgEPMx6gkbdvkRGRoZejo5qfRg5JSVFuo7W29sbCoUCgwcPfuUAiIiIyjqtk21+fr50Awng2WFbpVJpkKCIiIjKEq3P2QohEBISIp0jffLkCUJDQ9US7tatW/UbIRERUSmndbINDg5Wme7bt6/egyEiIiqLtE62q1evNmQcREREZVaxHh5PRERE2mOyJSIiMjAmWyIiIgNjsiUiIjIwJlsiIiIDY7IlIiIyMCZbIiIiA2OyJSIiMjAmWyIiIgNjsiUiIjIwJlsiIiIDY7IlIiIyMCZbIiIiA2OyJSIiMjAmWyIiIgNjsiUiIjIwJlsiIiIDY7IlIiIyMCZbIiIiAyt1yTY7Oxv16tWDTCbD6dOnVeYlJycjMDAQSqUS9vb2GDlyJHJyckomUCIiov/PpKQD0NWECRPg6uqKM2fOqJTn5eWhU6dOcHBwwOHDh5Geno7g4GAIIbBo0aISipaIiKiU7dnu2bMH0dHRmDNnjtq86OhonD9/Hj/88APq16+P1q1bY+7cuVi+fDkyMzNLIFoiIqJnSk2yvXnzJoYMGYLvv/8eFhYWavPj4uLg4+MDV1dXqaxdu3bIzs7GiRMnCm03OzsbmZmZKi8iIiJ9KhXJVgiBkJAQhIaGws/PT2OdtLQ0ODk5qZTZ2trCzMwMaWlphbYdFRUFGxsb6eXm5qbX2ImIiEo02UZEREAmkxX5On78OBYtWoTMzEyEh4cX2Z5MJlMrE0JoLC8QHh6OjIwM6XXt2rVX3i4iIqLnlegAqeHDh6N3795F1vH09MTUqVNx9OhRyOVylXl+fn7o06cP1q5dC2dnZxw7dkxl/r179/D06VO1Pd7nyeVytXaJiIj0qUSTrb29Pezt7V9ab+HChZg6dao0nZKSgnbt2mHTpk1o1KgRAMDf3x/Tpk1DamoqXFxcADwbNCWXy+Hr62uYDSAiItJCqbj0x93dXWXa0tISAFCpUiVUrFgRANC2bVvUrFkT/fr1w+zZs3H37l2MGzcOQ4YMgbW19WuPmYiIqECpGCClDWNjY+zatQsKhQJNmzZFz5490bVrV42XCREREb1OpWLP9kWenp4QQqiVu7u7Y+fOnSUQERERUeHKzJ4tERHRm4rJloiIyMCYbImIiAyMyZaIiMjAmGyJiIgMjMmWiIjIwJhsiYiIDIzJloiIyMCYbImIiAyMyZaIiMjAmGyJiIgMjMmWiIjIwJhsiYiIDIzJloiIyMCYbImIiAyMyZaIiMjASuXD4w2p4KH0Dx89KeFIiIiopBTkgIKc8KpkQl8tlRFXrlxBpUqVSjoMIiJ6A1y+fBne3t6v3A73bF9gZ2cHAEhOToaNjU0JR1M6ZGZmws3NDdeuXYO1tXVJh1NqsN90xz4rHvab7jIyMuDu7i7lhFfFZPsCI6Nnp7FtbGz4ptSRtbU1+6wY2G+6Y58VD/tNdwU54ZXb0UsrREREVCgmWyIiIgNjsn2BXC7H5MmTIZfLSzqUUoN9VjzsN92xz4qH/aY7ffcZRyMTEREZGPdsiYiIDIzJloiIyMCYbImIiAyMyZaIiMjAmGwBJCUlYdCgQfDy8oK5uTkqVaqEyZMnIycnR6VecnIyAgMDoVQqYW9vj5EjR6rVedtMmzYNTZo0gYWFBcqVK6exDvtN3eLFi+Hl5QWFQgFfX18cOnSopEN6oxw8eBCBgYFwdXWFTCbD9u3bVeYLIRAREQFXV1eYm5ujRYsWOHfuXMkE+4aIiopCw4YNYWVlBUdHR3Tt2hUXL15UqcN+U/Xdd9+hTp060s0+/P39sWfPHmm+PvuLyRbAhQsXkJ+fj6VLl+LcuXOYP38+lixZgs8++0yqk5eXh06dOiErKwuHDx/Gxo0b8dNPP+GTTz4pwchLXk5ODnr06IFhw4ZpnM9+U7dp0yaMHj0an3/+OU6dOoX33nsPHTp0QHJyckmH9sbIyspC3bp18c0332icP2vWLMybNw/ffPMN4uPj4ezsjDZt2uDBgwevOdI3R2xsLD7++GMcPXoUMTExyM3NRdu2bZGVlSXVYb+pqlixImbMmIHjx4/j+PHjaNmyJbp06SIlVL32lyCNZs2aJby8vKTp3bt3CyMjI3Hjxg2pbMOGDUIul4uMjIySCPGNsnr1amFjY6NWzn5T984774jQ0FCVsurVq4tPP/20hCJ6swEQ27Ztk6bz8/OFs7OzmDFjhlT25MkTYWNjI5YsWVICEb6Zbt26JQCI2NhYIQT7TVu2trZixYoVeu8v7tkWIiMjQ+UG1HFxcfDx8YGrq6tU1q5dO2RnZ+PEiRMlEWKpwH5TlZOTgxMnTqBt27Yq5W3btsWRI0dKKKrSJTExEWlpaSp9KJfL0bx5c/bhczIyMgD893AV9lvR8vLysHHjRmRlZcHf31/v/cVkq8Hly5exaNEihIaGSmVpaWlwcnJSqWdrawszMzOkpaW97hBLDfabqjt37iAvL0+tT5ycnN7K/iiOgn5iHxZOCIGxY8fi3XffhY+PDwD2W2HOnj0LS0tLyOVyhIaGYtu2bahZs6be+6tMJ9uIiAjIZLIiX8ePH1dZJiUlBe3bt0ePHj0wePBglXkymUxtHUIIjeWlWXH6rShvS7/p4sVtf9v7ozjYh4UbPnw4EhISsGHDBrV57DdV1apVw+nTp3H06FEMGzYMwcHBOH/+vDRfX/1Vph+xN3z4cPTu3bvIOp6entLfKSkpCAgIgL+/P5YtW6ZSz9nZGceOHVMpu3fvHp4+far2y6e007XfivI29Zs27O3tYWxsrPbL+NatW29lfxSHs7MzgGd7ai4uLlI5+/CZESNG4Oeff8bBgwdRsWJFqZz9ppmZmRkqV64MAPDz80N8fDy+/vprTJw4EYD++qtM79na29ujevXqRb4UCgUA4MaNG2jRogUaNGiA1atXqz3D0N/fH3/99RdSU1OlsujoaMjlcvj6+r7W7TI0XfrtZd6mftOGmZkZfH19ERMTo1IeExODJk2alFBUpYuXlxecnZ1V+jAnJwexsbFvdR8KITB8+HBs3boV+/btg5eXl8p89pt2hBDIzs7Wf3/pYfBWqXfjxg1RuXJl0bJlS3H9+nWRmpoqvQrk5uYKHx8f0apVK3Hy5Enx22+/iYoVK4rhw4eXYOQl7+rVq+LUqVMiMjJSWFpailOnTolTp06JBw8eCCHYb5ps3LhRmJqaipUrV4rz58+L0aNHC6VSKZKSkko6tDfGgwcPpPcSADFv3jxx6tQpcfXqVSGEEDNmzBA2NjZi69at4uzZs+LDDz8ULi4uIjMzs4QjLznDhg0TNjY24sCBAyrfYY8ePZLqsN9UhYeHi4MHD4rExESRkJAgPvvsM2FkZCSio6OFEPrtLyZb8eyyFQAaX8+7evWq6NSpkzA3Nxd2dnZi+PDh4smTJyUU9ZshODhYY7/t379fqsN+U/ftt98KDw8PYWZmJho0aCBdnkHP7N+/X+P7Kjg4WAjx7DKWyZMnC2dnZyGXy0WzZs3E2bNnSzboElbYd9jq1aulOuw3VQMHDpQ+hw4ODqJVq1ZSohVCv/3FR+wREREZWJk+Z0tERPQmYLIlIiIyMCZbIiIiA2OyJSIiMjAmWyIiIgNjsiUiIjIwJlsiIiIDY7IlIiIyMCZbojeQTCbD9u3bSzoMvQsJCZGeHPWq2/f806kWLFigl/iIDIXJlug1eT7RmJqawsnJCW3atMGqVauQn5+vUjc1NRUdOnTQqt3Slpjbt2+v0/YVZty4cUhNTVV5sg3Rm4rJlug1Kkg0SUlJ2LNnDwICAjBq1Ch07twZubm5Uj1nZ2fI5fISjNRw5HK5XrbP0tISzs7OMDY21lNkRIbDZEv0GhUkmgoVKqBBgwb47LPPsGPHDuzZswdr1qyR6j2/t5qTk4Phw4fDxcUFCoUCnp6eiIqKAvDfc4W7desGmUwmTV++fBldunSBk5MTLC0t0bBhQ/z2228qsXh6emL69OkYOHAgrKys4O7urvYc5+vXr6N3796ws7ODUqmEn5+fyvOJf/nlF/j6+kKhUMDb2xuRkZEqPxq0kZSUBJlMhs2bN+O9996Dubk5GjZsiH/++Qfx8fHw8/ODpaUl2rdvj9u3b+vUNtGbgsmWqIS1bNkSdevWxdatWzXOX7hwIX7++Wds3rwZFy9exA8//CAl1fj4eADA6tWrkZqaKk0/fPgQHTt2xG+//YZTp06hXbt2CAwMRHJyskrbc+fOhZ+fH06dOoWwsDAMGzYMFy5ckNpo3rw5UlJS8PPPP+PMmTOYMGGCdMj7119/Rd++fTFy5EicP38eS5cuxZo1azBt2rRi9cPkyZPxxRdf4OTJkzAxMcGHH36ICRMm4Ouvv8ahQ4dw+fJlTJo0qVhtE5U4/TyoiIheJjg4WHTp0kXjvF69eokaNWpI0wDEtm3bhBBCjBgxQrRs2VLk5+drXPb5ukWpWbOmWLRokTTt4eEh+vbtK03n5+cLR0dH8d133wkhhFi6dKmwsrIS6enpGtt77733xPTp01XKvv/+e+Hi4lJoDJr6IDExUQAQK1askMo2bNggAIjff/9dKouKihLVqlVTa9PDw0PMnz+/0HUSvQlMSjbVExEACCEgk8k0zgsJCUGbNm1QrVo1tG/fHp07d0bbtm2LbC8rKwuRkZHYuXMnUlJSkJubi8ePH6vt2dapU0f6WyaTwdnZGbdu3QIAnD59GvXr14ednZ3GdZw4cQLx8fEqe7J5eXl48uQJHj16BAsLC622XVMsTk5OAIDatWurlBXERlTaMNkSvQH+/vtveHl5aZzXoEEDJCYmYs+ePfjtt9/Qs2dPtG7dGlu2bCm0vfHjx+PXX3/FnDlzULlyZZibm6N79+7IyclRqWdqaqoyLZPJpMPE5ubmRcacn5+PyMhIfPDBB2rzFApFkctq8nwsBT88Xix7cdQ2UWnBZEtUwvbt24ezZ89izJgxhdaxtrZGr1690KtXL3Tv3h3t27fH3bt3YWdnB1NTU+Tl5anUP3ToEEJCQtCtWzcAz86/JiUl6RRXnTp1sGLFCmk9L2rQoAEuXryIypUr69Qu0duIyZboNcrOzkZaWhry8vJw8+ZN7N27F1FRUejcuTP69++vcZn58+fDxcUF9erVg5GREX788Uc4OzujXLlyAJ6NKv7999/RtGlTyOVy2NraonLlyti6dSsCAwMhk8nw5Zdf6rxX+OGHH2L69Ono2rUroqKi4OLiglOnTsHV1RX+/v6YNGkSOnfuDDc3N/To0QNGRkZISEjA2bNnMXXq1FftKqIyhaORiV6jvXv3wsXFBZ6enmjfvj3279+PhQsXYseOHYVeL2ppaYmZM2fCz88PDRs2RFJSEnbv3g0jo2cf37lz5yImJgZubm6oX78+gGcJ2tbWFk2aNEFgYCDatWuHBg0a6BSrmZkZoqOj4ejoiI4dO6J27dqYMWOGFGe7du2wc+dOxMTEoGHDhmjcuDHmzZsHDw+PV+ghorJJJoQQJR0EEb0dQkJCcP/+fb3e8crT0xOjR4/G6NGj9dYmkb5xz5aIXqudO3fC0tISO3fufKV2pk+fDktLS7UR1kRvIu7ZEtFrc+vWLWRmZgIAXFxcoFQqi93W3bt3cffuXQCAg4MDbGxs9BIjkSEw2RIRERkYDyMTEREZGJMtERGRgTHZEhERGRiTLRERkYEx2RIRERkYky0REZGBMdkSEREZGJMtERGRgf0/RJ0rgH8XgZIAAAAASUVORK5CYII=", "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, 0), 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, -47.87),\n", " width=50,\n", " height=47.87,\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, -(16.77 + 1.52)),\n", " width=2,\n", " height=(16.77 + 1.52),\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, 0),\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, -(16.77 + 1.52)),\n", " width=2,\n", " height=1.52,\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=3, color=\"#00035b\")\n", "ax.add_patch(pumping_arrow)\n", "ax.text(x=0.5, y=11, s=r\"$ D = 0.671$ m\", fontsize=\"large\")\n", "\n", "# last line\n", "line = plt.Line2D(xdata=[-200, 1200], ydata=[0, 0], color=\"k\")\n", "ax.add_line(line)\n", "\n", "\n", "ax.set_xlim([-20, 30])\n", "ax.set_ylim([-47, 20])\n", "ax.set_xlabel(\"Distance [m]\")\n", "ax.set_ylabel(\"Relative height [m]\")\n", "ax.set_title(\"Conceptual Model - Pratt County Example\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 2. Set basic parameters" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "rw = 0.125 # well radius, m\n", "rc = 0.064 # well casing radius, m\n", "L = 1.52 # screen length, m\n", "b = -47.87 # aquifer thickness, m\n", "zt = -16.77 # depth to top of screen, m\n", "H0 = 0.671 # initial displacement in the well, m\n", "zb = zt - L # bottom of screen, m" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 3. Calculate the added volume\n", "\n", "As we will see later, the input for the slug test in TTim is the added or removed volume. Therefore we must first convert our measured displacement into volume." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "slug: 0.00863 m^3\n" ] } ], "source": [ "Q = np.pi * rc**2 * H0\n", "print(\"slug:\", round(Q, 5), \"m^3\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 4. Load data from well" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Data from the slug well is loaded from a text file, where the first column is the time in seconds and the second column is the head." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "data = np.loadtxt(\"data/slug.txt\", skiprows=1)\n", "to = data[:, 0] / 60 / 60 / 24 # convert time to days\n", "ho = data[:, 1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 5. Create a conceptual model in TTim\n", "\n", "We conceptualize the aquifer as a three-layer model, one layer above the screen, one layer at the screen top and bottom and another layer just below it.\n", "\n", "We use ```Model3D``` method to build this model. Details on how to set the model can be seen in notebook: [Unconfined - 1 - Vennebulten](unconfined1_vennebulten.ipynb).\n", "\n", "The setting of the slug well is slightly different from the pumping well. We detail the differences below:\n", "* the ```tsandQ``` argument in the ```Well``` object has a different meaning. Instead of meaning the time of start or shutdown and the pumping rate of the pumping well, it means the time a volume is instantaneously added or removed from the well. In our case, we defined it as ```[(0, -Q)]``` where ```0``` is the moment in time when we added the slug and ```Q``` is the added volume. A negative sign means a volume is added. Otherwise, it would mean an extracted volume.\n", "* the ```wbstype``` argument is set to ```' slug'```, so TTim knows the ```tsandQ``` argument means time and instant volumes, instead of pumping rates." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "self.neq 1\n", "solution complete\n" ] } ], "source": [ "ml = ttm.Model3D(kaq=10, z=[0, zt, zb, b], Saq=1e-4, kzoverkh=1, tmin=1e-6, tmax=0.01)\n", "w = ttm.Well(ml, xw=0, yw=0, rw=rw, rc=rc, tsandQ=[(0, -Q)], layers=1, wbstype=\"slug\")\n", "ml.solve()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 6. Model calibration\n", "\n", "We calibrate both hydraulic conductivity and specific storage, considering uniform parameters in all layers. Since the observations are inside the pumping well, we use the `seriesinwell` function to add the observations. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "............................\n", "Fit succeeded.\n" ] } ], "source": [ "ca = ttm.Calibrate(ml)\n", "ca.set_parameter(name=\"kaq0_2\", initial=10)\n", "ca.set_parameter(name=\"Saq0_2\", initial=1e-4)\n", "ca.seriesinwell(name=\"obs\", element=w, t=to, h=ho)\n", "ca.fit()" ] }, { "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_26.0490850.0231380.382507-infinf10.0000[6.049085184433308, 6.049085184433308, 6.04908...
Saq0_20.0002130.0000104.490863-infinf0.0001[0.00021324705409226575, 0.0002132470540922657...
\n", "
" ], "text/plain": [ " optimal std perc_std pmin pmax initial \\\n", "kaq0_2 6.049085 0.023138 0.382507 -inf inf 10.0000 \n", "Saq0_2 0.000213 0.000010 4.490863 -inf inf 0.0001 \n", "\n", " parray \n", "kaq0_2 [6.049085184433308, 6.049085184433308, 6.04908... \n", "Saq0_2 [0.00021324705409226575, 0.0002132470540922657... " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "display(ca.parameters)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we plot the results in heads over initial heads, as done for the graphical solution." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcoAAAE/CAYAAADGw4N2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABV5klEQVR4nO3dd3yN1x/A8c/N3kvWRZZNzUZoqBEqROlApdWKXUoRqzX6o1RRRWlrtEV1qBqlpU0R1KZitfZMBE1EjOxxkzy/P9LcujIkJLkZ3/frdV/ynHue537vPVe+Oec5z3lUiqIoCCGEECJPBvoOQAghhCjLJFEKIYQQBZBEKYQQQhRAEqUQQghRAEmUQgghRAEkUQohhBAFkEQphBBCFEASpRBCCFEASZRCCCFEASRRikJZtWoVKpUKlUrF7t27cz2vKAq1atVCpVLRvn37Yn1tlUrF+++/X+T9IiIiUKlUrFq1qlD1ch4GBgbY29vTsWNHtm/f/nhBF7P27dvrfK7Jycm8//77ebZFSevfv7/O55Xfo3///gB4enrSrVu3Uo/zcTzud01UbEb6DkCUL9bW1qxYsSJXMtyzZw9XrlzB2tpaP4EVg5EjR9KnTx8yMzM5f/4806dPp2vXruzatYu2bdvqOzwdycnJTJ8+HaDY/zB5lP/9738MGzZMu338+HFGjBjBrFmz8PPz05Y7OTmValxClBRJlKJIAgMDWb16NYsXL8bGxkZbvmLFCnx9fYmPj9djdE/G3d2dZ555BoDWrVtTu3Zt2rVrx4oVK8pcotSnmjVrUrNmTe12amoqALVr19Z+fsUlJSUFMzMzVCpVsR63PFIUhdTUVMzNzfUdSqUjQ6+iSF577TUA1qxZoy2Li4vjp59+YuDAgXnuc/fuXYYPH061atUwMTGhRo0aTJkyhbS0NJ168fHxDBkyhCpVqmBlZUWXLl24ePFinse8dOkSffr0wdnZGVNTU+rXr8/ixYuL6V1ma968OQC3bt3SKY+Ojmbo0KFUr14dExMTvLy8mD59OhkZGTr1li5dSpMmTbCyssLa2pp69eoxefJk7fPvv/9+ngkgZ5g7IiIiz7giIiK0vbXp06fnGuq8ffs2b775Jm5ubpiamuLk5ETr1q3ZsWPH434UxWLr1q08/fTTmJubU69ePVauXKnzfM773r59OwMHDsTJyQkLCwvt92Tt2rX4+vpiaWmJlZUVnTt35sSJE7le5+jRo7zwwgs4ODhgZmZGs2bNWLdu3WPFfPv2bYYPH06DBg2wsrLC2dmZDh06sG/fPm0dRVGoXbs2nTt3zrV/YmIitra2jBgxQlsWHx/P+PHj8fLywsTEhGrVqhEcHExSUpLOviqVirfffptly5ZRv359TE1N+eabb4BHf7dE8ZIepSgSGxsbevXqxcqVKxk6dCiQnTQNDAwIDAxk4cKFOvVTU1Px8/PjypUrTJ8+ncaNG7Nv3z5mz57NyZMn+e2334DsXzYvvfQSBw8eZOrUqfj4+HDgwAECAgJyxXD27FlatWqFu7s78+fPx9XVlW3btjFq1ChiY2OZNm1asbzX8PBwAOrUqaMti46OpkWLFhgYGDB16lRq1qzJoUOHmDlzJhEREXz99dcA/PjjjwwfPpyRI0cyb948DAwMuHz5MmfPnn3iuNRqNVu3bqVLly4MGjSIwYMHA/8Ndfbt25fjx4/z4YcfUqdOHe7fv8/x48e5c+fOE7/24/rrr78YN24cEydOxMXFheXLlzNo0CBq1aqVq7c+cOBAnn/+eb777juSkpIwNjZm1qxZvPfeewwYMID33nuP9PR0Pv74Y9q0acORI0do0KABAH/88QddunShZcuWLFu2DFtbW3788UcCAwNJTk7W/jFRWHfv3gVg2rRpuLq6kpiYyKZNm2jfvj07d+6kffv2qFQqRo4cSXBwMJcuXaJ27dra/b/99lvi4+O1iTI5OZl27dpx48YNJk+eTOPGjTlz5gxTp07l1KlT7NixQ+ePp59//pl9+/YxdepUXF1dcXZ2LtHvlsiHIkQhfP311wqghIWFKX/88YcCKKdPn1YURVF8fHyU/v37K4qiKE899ZTSrl077X7Lli1TAGXdunU6x/voo48UQNm+fbuiKIry+++/K4CyaNEinXoffvihAijTpk3TlnXu3FmpXr26EhcXp1P37bffVszMzJS7d+8qiqIo4eHhCqB8/fXXBb63nHofffSRotFolNTUVOXkyZOKr6+volarlfDwcG3doUOHKlZWVsq1a9d0jjFv3jwFUM6cOaONxc7OrsDXnTZtmpLXf8Gcz/rB123Xrp3O53r79u1cn0sOKysrJTg4uMDXLk4534f169fn+byHh4diZmam85mlpKQoDg4OytChQ7VlOe87KChIZ//IyEjFyMhIGTlypE55QkKC4urqqvTu3VtbVq9ePaVZs2aKRqPRqdutWzdFrVYrmZmZBb6X/D7THBkZGYpGo1E6duyovPzyy9ry+Ph4xdraWhk9erRO/QYNGih+fn7a7dmzZysGBgZKWFiYTr0NGzYogBISEqITi62trfb7nKMw3y1RvGToVRRZu3btqFmzJitXruTUqVOEhYXlO+y6a9cuLC0t6dWrl055zl/2O3fuBLJ7AgCvv/66Tr0+ffrobKemprJz505efvllLCwsyMjI0D66du1Kamoqhw8ffqz39e6772JsbIyZmRlNmzbl9OnTbNmyBU9PT22dX3/9FT8/P6pWrarz2jk93z179gDQokUL7t+/z2uvvcYvv/xCbGzsY8X0OFq0aMGqVauYOXMmhw8fRqPRFGq/B99PRkYGSjHeqrZp06a4u7trt83MzKhTpw7Xrl3LVbdnz54629u2bSMjI4OgoCCd+MzMzGjXrp125u/ly5c5f/689jv08HcjKiqKCxcuFDn2ZcuW8fTTT2NmZoaRkRHGxsbs3LmTc+fOaetYW1szYMAAVq1apR1C3bVrF2fPnuXtt9/W1vv1119p2LAhTZs21Ymvc+fOec4o79ChA/b29jpl+vxuVVaSKEWRqVQqBgwYwPfff8+yZcuoU6cObdq0ybPunTt3cHV1zXUuztnZGSMjI+1w4J07dzAyMqJKlSo69VxdXXMdLyMjg88++wxjY2OdR9euXQEe+xfH6NGjCQsLY//+/cybNw+NRsOLL76oM2R569YttmzZkuu1n3rqKZ3X7tu3LytXruTatWv07NkTZ2dnWrZsSWho6GPFVhRr166lX79+LF++HF9fXxwcHAgKCiI6OjrffSIiInK9p5ykXxweblcAU1NTUlJScpWr1Wqd7ZxzxD4+PrliXLt2rfYzz6k3fvz4XPWGDx8OFP27sWDBAt566y1atmzJTz/9xOHDhwkLC6NLly65Yh85ciQJCQmsXr0agM8//5zq1avz4osv6ryXv//+O1d81tbWKIqSK76HPwvQ73erspJzlOKx9O/fn6lTp7Js2TI+/PDDfOtVqVKFP//8E0VRdJJlTEwMGRkZODo6autlZGRw584dnV+qD/9yt7e3x9DQkL59++pMkHiQl5fXY72n6tWrayfwtG7dGldXV9544w2mTZvG559/DoCjoyONGzfO9z1XrVpV+/OAAQMYMGAASUlJ7N27l2nTptGtWzcuXryIh4cHZmZmAKSlpWFqaqrd70l7CI6OjixcuJCFCxcSGRnJ5s2bmThxIjExMWzdujXfuMPCwnTK6tat+0RxPK6H/6jK+Y5s2LABDw+PfPfLqTdp0iR69OiRZ52ivqfvv/+e9u3bs3TpUp3yhISEXHVr1apFQEAAixcvJiAggM2bNzN9+nQMDQ11YjQ3N881kenh95Ajv9m+j/puieIliVI8lmrVqjFhwgTOnz9Pv3798q3XsWNH1q1bx88//8zLL7+sLf/222+1zwP4+fkxd+5cVq9ezahRo7T1fvjhB53jWVhY4Ofnx4kTJ2jcuDEmJibF+bZ0vP766yxfvpyvvvqKCRMm4OHhQbdu3QgJCaFmzZq5hsTyY2lpSUBAAOnp6bz00kucOXMGDw8P7ZDu33//jY+Pj7b+li1bHnnMnMSaV4/sQe7u7rz99tvs3LmTAwcO5FvPxMRE+0dCWdO5c2eMjIy4cuVKrmHZB9WtW5fatWvz119/MWvWrGJ5bZVKpfNHDGS316FDh3Bzc8tVf/To0fj7+9OvXz8MDQ0ZMmSIzvPdunVj1qxZVKlS5bH/oHtQft8tUbwkUYrHNmfOnEfWCQoKYvHixfTr14+IiAgaNWrE/v37mTVrFl27duW5554DwN/fn7Zt2/LOO++QlJRE8+bNOXDgAN99912uYy5atIhnn32WNm3a8NZbb+Hp6UlCQgKXL19my5Yt7Nq1q9je40cffUTLli354IMPWL58OTNmzCA0NJRWrVoxatQo6tatS2pqKhEREYSEhLBs2TKqV6/OkCFDMDc3p3Xr1qjVaqKjo5k9eza2trbapNi1a1ccHBwYNGgQM2bMwMjIiFWrVnH9+vVHxmVtbY2Hhwe//PILHTt2xMHBAUdHR+zt7fHz86NPnz7Uq1cPa2trwsLC2Lp1a769rLLO09OTGTNmMGXKFK5evUqXLl2wt7fn1q1bHDlyBEtLS+3iC1988QUBAQF07tyZ/v37U61aNe7evcu5c+c4fvw469evL9Jrd+vWjQ8++IBp06bRrl07Lly4wIwZM/Dy8sp1ORBAp06daNCgAX/88QdvvPEGzs7OOs8HBwfz008/0bZtW8aMGUPjxo3JysoiMjKS7du3M27cOFq2bFlgTIX5bolipufJRKKceHDWa0EenvWqKIpy584dZdiwYYparVaMjIwUDw8PZdKkSUpqaqpOvfv37ysDBw5U7OzsFAsLC6VTp07K+fPn85yJGB4ergwcOFCpVq2aYmxsrDg5OSmtWrVSZs6cqVOHIsx6/fjjj/N8/pVXXlGMjIyUy5cvK4qSPeN01KhRipeXl2JsbKw4ODgo3t7eypQpU5TExERFURTlm2++Ufz8/BQXFxfFxMREqVq1qtK7d2/l77//1jn2kSNHlFatWimWlpZKtWrVlGnTpinLly9/5KxXRVGUHTt2KM2aNVNMTU0VQOnXr5+SmpqqDBs2TGncuLFiY2OjmJubK3Xr1lWmTZumJCUlFfg5PK7CzHp9/vnnc5U//J4e9R37+eefFT8/P8XGxkYxNTVVPDw8lF69eik7duzQqffXX38pvXv3VpydnRVjY2PF1dVV6dChg7Js2bJHvpeHv2tpaWnK+PHjlWrVqilmZmbK008/rfz8889Kv379FA8PjzyP8f777yuAcvjw4TyfT0xMVN577z2lbt26iomJiWJra6s0atRIGTNmjBIdHa0Ty4gRI3LtX9jvlig+KkUpxqltQghRyTVv3hyVSpXrnK8ov2ToVQghnlB8fDynT5/m119/5dixY2zatEnfIYliJIlSCCGe0PHjx/Hz86NKlSpMmzaNl156Sd8hiWIkQ69CCCFEAfS64MDevXvp3r07VatWRaVS8fPPPz9ynz179uDt7Y2ZmRk1atRg2bJlJR+oEEKISkuviTIpKYkmTZpoL+Z+lPDwcLp27UqbNm04ceIEkydPZtSoUfz0008lHKkQQojKqswMvapUKjZt2lTg2P67777L5s2bddZYHDZsGH/99ReHDh0qhSiFEEJUNuVqMs+hQ4fw9/fXKevcuTMrVqxAo9FgbGyca5+0tDSd+x5mZWVx9+5dqlSpIjeDFUKISkxRFBISEqhatSoGBvkPsJarRBkdHY2Li4tOmYuLCxkZGcTGxua5gPDs2bO1q3YIIYQQD7t+/TrVq1fP9/lylSgh9yLBOSPH+fUOJ02axNixY7XbcXFxuLu7Ex4ejrW19SNfT6PR8Mcff+Dn56ftsf4SFkGvgOJdF1MxNART0+yHuTmYmYGJCUrOz2Zm/JMGp2LTSDEyIc3YhGZ1q1Lbwym7voVFdt1/f8bcHKysssssLf97mJlBOe1J59UWQj+kLcoOaYvHl5CQgJeX1yNzQblKlK6urrnuJhETE5Pn7ZlymJqa5lrUGMDBwQEbG5tHvqZGo8HCwoIqVapgbGxMVFwKc3dex6RpAEZZmRhnajDNyuS5mnaYkQVpaZCe/t+/OT/n9cjK+u+FMjMhOTn7ce9enrHYAQ0eLPjrkeHnZmhIlqUVGnMLDG1tMLKxBhsbsLbOftjYZD9sbf/719aWWCNzbmQZo/ZQ4+Khzk66DyXcqLgUwmOT8HK0RG1rnuulH/X8ozzcFkJ/pC3KDmmLx5fzeT3qNFy5SpS+vr657qywfft2mjdvXmpfkPDYJFINjXmvs+4tntYMeQbfmnkn6zwpCmRk/Jc0U1Oz/01Jyf45NfW/n1NSuHQthq93nMNMk45ZRhqmGemYa9J4oY49rsZZkJSUnWQf/DcpCRITs//NuctEZiYG8XGYxsfBrahCh+v470PL0BDs7bWPKENzjsTBXXMbwsysaPNMXZ72rgNVqoCjI1tupvHOnihSjEwxUMHsHo0I9PnvRr5PmkSLQ1mIQQhR9ug1USYmJnL58mXtdnh4OCdPnsTBwQF3d3cmTZrEzZs3tbdkGjZsGJ9//jljx45lyJAhHDp0iBUrVrBmzZpSi9nL0RIDFWQ9MFfYUKXC09GiaAdSqcDYOPthZfXI6lZxKfx4Z1eu1+0+0Q8K80s9M5PoqDu8+NE2LNJSsEhPwSo9BWtNKh8H1MA+IxUSEiAuDuLjsx9xcaTduceFCzewTkvEOi0Z29REjLMys3vAsbHZD0ANvPjg6z10R6fu/z6SjM24Y2HLnW/tSG1SC7Nqas5kmrHuhoYYC3vuWNkT1NOXbp29s4eQS8nasEgmbTxFlkKeiVwIUXnpNVEePXoUPz8/7XbOucR+/fqxatUqoqKiiIyM1D7v5eVFSEgIY8aMYfHixVStWpVPP/20wHvUFTe1rTmzezRi8sbTZCoKhioVs3o0LPEeyBO/rqEhV9MMuGXpAJa6T51vl39v+NiVWPp89ed/BYqCWUYaq3vUxdsWuHePc2evsfKXo9ikJmKfmoB9Sjx2KQk8a6/CNimO9FsxEHsHk6wMLDWpWMal4h53C/65AMBTgM50q9X//mtri8ZVTZKjKyYe1TF1r45HXBz3kjX8Y+uM61O1cfXIPYGrqD3DqLgUbZKE7D+CJv10inqu1jRxK9w9J4UQFZdeE2X79u0p6DLOVatW5Spr164dx48fL8GoHi3Qx522dZyIiE3G09Gi1IbpnvR1H6c3nGsflQqNiTlVG9bU9mTtGjfnpwj7XMfdP9EPW1tz7sSl0Hr2TizTknFIjqNKUhzOKXHMaetCfMRN/th7CqfEezgl3cc56S7Oifcwz0iDuDiM4+Kwu3Be20NtCrB0KTl3+Uu3tMakhid4eICnJycM7fgiMotrdq5ct1fzv1db0LaOU4GJMzw2SSd2gCzgpSUHmSM9S1HGZGVlkZ6ert3WaDQYGRmRmppKZmamHiMre4yNjTE0NHzi45Src5RlidrWXC/nsZ7kdR+nV1qYfR5VR21rzuyejZm88TQJppbccKjGrB4NsfNxJyUuhelzHhpSBn7u25CxC3/HKeEOrol3cE24gzohNvsRn/2vQ0o8JkkJcOpU9gNoBjy4qGHsMjuu2bsSZafmkEM1fLs8Q+1nn+aKfVU83J1R25rn+QcEZJ9GnrzxNG3rOKG2NZdzmELv0tPTCQ8PJ+uBiYCKouDq6sr169fl2vA82NnZ4erq+kSfjSTKSuZxeqWF2edRdfJ7Pr8km2BmwaUqblyq4pZvXObpqVSNv83nrRyon36XGyfOcnLvSarH3cLtfjRVUuJxTL6PY/J9vG+ez95p3/cAOAE3bJyJalAf9TPNWGvhykfXDLno6Ea82X/njDMVhYjYZPZevJ3rHOajeqpCFCdFUYiKisLQ0BA3NzftBfJZWVkkJiZiZWVV4EXzlY2iKCQnJxMTEwOQ53X2hSWJshJ6nF5pYfZ5VJ38ns8riUbFpeTq5WX/Paig/PtTiokZEU7u2PXMntBkGJfCqAd6p9ZpSbjfi8LjfjSe9/7B6+4/eN27SY27N3FIiad6fAwcjoHDe/ABNvz7Ojetnbjg5MF5Z08uONfAPtKRvltvkaXKHsLJUmDiT6dQ/RufTP4RpSEjI4Pk5GSqVq2KxQMT3XKGYs3MzCRRPsTcPPv3TUxMDM7Ozo89DCuJUpQJDyfRvHqaH7xYn7//PsW6cEOyFPIe4n1gnyRTS8661uKMa61cr2efHEeNuzepdec6o9Uaqv4TDmfOwI0bVEu4TbWE23S4ejS78ua5/G1syjknL0651uKUa21OudbkchU3MMiO5cEhWiFKQs75RxMTEz1HUr7k/FGh0WgkUYqK5+GepqOFEZa3/mZ4j7bcjEsv1BDv3ou3tYnTAFD+fdyzsOWYhS0n3Z4i+MFLbO7fJ/bwMeKPHMfl2kUsz54m6++/sUhOxvuf83j/c177WknGZpxyrcWJqvU4WbUON895on6mgZzLFCVKzkMWTXF8XpIoRZn2YE9To9H8W2aGu2P+S049uE9BiTPPyUx2djh26Yhjl47aIoPMTEI27mH7d1t5KvoSjaMv89StK1ilp/DM9dM8c/10dsVNs0io7sFBu5qEVWtAmHtD3hzchcAWHpI8hSjHJFGKCq+gxFmopGVoSNdXOtDM35eI2GTcHS0IORfNyuVbafTPebz/uUCXpGvYXrmA9Y1r9LxxjZ6ndwEQ84M9V1r48rWJJ/vdmxDpUJXZPRvL+Uwh/rV79278/Py4d+8ednZ2+g4nT5IoRaXzuJfYPLhf72e8aFO/vzbh2tma8+fJqyyZ9T3eN87S4sYZmv1zAeekezj/EcLMf49xw8aJg7835V5wX+xffB4cHIrxnQkhSoIkSiEe08MJ191Lzb6a3uyp4Q2AaUY6TaMu0jLyFK2u/UWzf85TPf42vf8OhYGhKAYGJDZ+GgK6YN3rZaJq1CP8TrIMzwpRxshcYiGKSc6sW8N/Jw9kGJvS4c1efP7sa7zaZw5NRv9I394zWN7iZe541UGVlYX1yaNYz54J3t6o3NyI6NmXSW9+zPqD2WsgR8WlcPBKLFFxKfp8a6KCKe3vVVpaGqNGjcLZ2RkzMzOeffZZwsLCdOocOHCAJk2aYGZmRsuWLTn17yIiANeuXaN79+7Y29tjaWnJU089RUhISKnEDtKjFKJY5XUO1M7CmMkbT5NqbMbBGt48O7wPPlvP4xwXS7vwY3S8EsazESdwTbxLn7+20uevrST88hGX2jzHAuuG7PLyRmNiKtdqimKhjxsAvPPOO/z000988803eHh4MHfuXDp37qxzU4wJEyawaNEiXF1dmTx5Mi+88AIXL17E2NiYESNGkJ6ezt69e7G0tOTs2bNYFeJmEsVFEqUQxezhIdmHk2fO2rLRNo6sbdKZtU06Y5qRju+1v3nu8p/4XzqMc9I9rHdsZimbSTQxJ7RWS0IutaXt0rGonWz1+O5EeZbXDQBK+hrgpKQkli5dyqpVqwgICADgq6++IjQ0lBUrVuDj4wPAtGnT6NSpEwDffPMN1atXZ9OmTfTu3ZvIyEh69uxJo0aNAKhRo0aJxJofSZRClIKHk+fDqw6lGZmwu2Zzdtdszv/838L75gX8Lx6k64X9VI+/zctnd/Py2d1odn1OUs9eXPF/CadO7VDbld6tyET5l9cNAHKWaSypRHnlyhU0Gg2tW7fWlhkbG9OiRQvOnTunTZS+vr7a5x0cHKhbty7nzp0DYNSoUbz11lts376d5557jp49e9K4ceMSiTcvco5SiFL28LlMQ5WKnk9X024bGBjSacjLzOk4iGeHraTHGx+z0vsFblk5YHzvLpbLv6Rx766k1qjNX29PhKgoOZcpCiXnBgAPeqz76RZBzh2iHr7wX1GURy4GkPP84MGDuXr1Kn379uXUqVM0b96czz77rGQCzoMkSiH0INDHnf0T/Vgz5Bn2T/Rjfu+mOttD29XMTqYGBhyvVp8POw3ll18OERT4ARuf8iPZ2BSve//QZPFHZLm5caZ5e5a/+xltZoWyNizy0QGISimvP9JK+n66tWrVwsTEhP3792vLNBoNR48epX79+tqyw4cPa3++d+8eFy9epF69etoyNzc3hg0bxsaNGxk3bhxfffVVicX8MBl6FUJP8lrf9lHnNmd5NmOvZzPeSx/O8+f30/vvUHxunuW5y0d47vIRbtg488OfAdxaPh2XWh76eFuijCvt++laWlry1ltvMWHCBBwcHHB3d2fu3LkkJyczaNAg/vrrLwBmzJhBlSpVcHFxYcqUKTg6OvLSSy8BEBwcTEBAAHXq1OHevXvs2rVLJ8mWNEmUQpRh+Z3bTDYxZ33jTvzUuBNed64T+Nd2Xjm1g+rxMbyz5xuyGqwhuVdvLvYZhEublnJdptBR2vfTnTNnDllZWfTt25eEhASaN2/Otm3bsLe316kzevRoLl26RJMmTdi8ebN2AfjMzExGjBjBjRs3sLGxoUuXLnzyySelFr9KyRlAriTi4+OxtbUlLi4OGxubR9bXaDSEhITQtWtXjI2NSyFCkR9pi+yp/Q+uVftOl7p8tPU8WQqYatLodn4/fU/8RtOoi9p9Dro3JmN0MG2D+0Ex3YZJ2qL0paamEh4ejpeXF2ZmZtryrKws4uPjsbGxkdts5SG/zw0Knw+kRylEOVLQdZppxqb83Pg56rwzgukrf2HgkZ8JuHCAVpF/w7iBaFbMJ2n0WM61DcBTbSe9TCEKSRKlEOVMYc5lzlbXZeSL76KOv03/Y1voc/J3rM+ewW7oINxsnPnc9xWaTg3mlVa15M4mQjyCJEohKoD8zmVG2Tgx228gS3x788bJEPqHbaZ6fAwfblvMzUPr2DtkFEMMGpJmYFxqq7QIUd7IgLYQFUxelwAEdmrE4md68+ywFbzf8U1uWTlQLf42bef/j51fvEnPUzshM5PJG0/LtZhCPER6lEJUQA8PxwIs3x9OmrEpq5q/wJomnXnt7+28dXg91eNvMz/kE9488hMft+1HxO0WADIcK8S/JFEKUUE9PBw7u0cj7YzZDBMz1O9NoP0Wf/oe+5Xhh9dTNzaS5Rs/4J8rWxna7HX+dqklw7FCIIlSiEoj3xmzxmb82KQzw45sZMjxzVT9K4zNf4XxU8MOzG0bVOKLZgtR1kmiFKISKXjG7Av89fdFIt8aQ48zf9Dz9C4CLhzg01avce2fZpIoRaUlk3mEqOTUtub41qyC2tacao3rML77OF7sO5+j1epjoUlj4p5VNH/JjzubQ2ThdVEpSaIUQmjlzJg9Xa0evV6fy/jnx5Lq4IjRxQtUefF5Yl7oRfepm1h/7Ia+QxUVTPv27QkODtZ3GHmSoVchhA7d4diO3L8/nq09h9L3+G+8dHYPbcNPMPPqYJq80YaouFRuxMXJ7FiRr/bt29O0aVMWLlwIwO7du/Hz8+PevXvY2dlp623cuLHMLocoiVIIkcuD5zIPxibx/nND+empDny09VMaxISz4NcFnL70B6/tvc9NayeZHSuemIODg75DyJcMvQohCpRzs99T6tq8EPQJH7XrR5qhMQ0vnOD35SPoeWonWVmKLFYgcunfvz979uxh0aJFqFQqVCoVfn5+ANjb26NSqejfvz+Qe+jV09OTmTNnEhQUhJWVFR4eHvzyyy/cvn2bF198ESsrKxo1asTRo0dL/H1IohRCFOjBlX4yDI340rc3yz/fyPGqdbFJT2Z+yCd8tXEm9ol3ORZxTyb8lBZFgaQk/TwKedOpRYsW4evry5AhQ4iKiuLGjRts2LABgAsXLhAVFcWiRYvy3f+TTz6hdevWnDhxgueff56+ffsSFBTEG2+8wfHjx6lVqxZBQUGU9E2wZOhVCPFID1+DqdFk0OHqRww5sokx+1bT6fKfNFt5nnduBbOrho8MxZaG5GSwssIAsCvt105MBEvLR1aztbXFxMQECwsLXF1dAahSpQoAzs7OOuco89K1a1eGDh0KwNSpU1m6dCk+Pj688sorALz77rv4+vpy69Yt7fFLgvQohRCF8uBlJGpbM3rVUvGF7yt077+Qc06eOCbHsXL9dKbu+BIjjUaGYsUTa9y4sfZnFxcXABo1apSrLCYmpkTjkB6lEOKx+LooDO/RlptxzxA+LIDDbwUz4NgWBh7bzDPXTzHixYlExCbLbNiSYmEBiYn6uXGzhUWpvMyDs2BV/y7yn1dZVlZWicah9x7lkiVLtHee9vb2Zt++fQXWX716NU2aNMHCwgK1Ws2AAQO4c+dOKUUrhHiQ2tYM35pVaFZHzQedhjKg1zTumNvQICaczd8EU+/ANn2HWHGpVNnDn/p4/JugCsPExITMzEydbUCnrKzTa6Jcu3YtwcHBTJkyhRMnTtCmTRsCAgKIjIzMs/7+/fsJCgpi0KBBnDlzhvXr1xMWFsbgwYNLOXIhxINyJvzsrdWCLgM/50+3hlinp2Df73UYO5ao2HiZ5FNJeXp68ueffxIREUFsbCweHh6oVCp+/fVXbt++TWJior5DfCS9JsoFCxYwaNAgBg8eTP369Vm4cCFubm4sXbo0z/qHDx/G09OTUaNG4eXlxbPPPsvQoUNLZXqwEKJggT7u7J/ox6djuuJ+4iC88072E598ws1mvoxauJXWc3axNizvP4RFxTR+/HgMDQ1p0KABTk5OaDQapk+fzsSJE3FxceHtt9/Wd4iPpLdzlOnp6Rw7doyJEyfqlPv7+3Pw4ME892nVqhVTpkwhJCSEgIAAYmJi2LBhA88//3y+r5OWlkZaWpp2Oz4+HgCNRoNGo3lknDl1ClNXlCxpi7Ijv7ZwtDDC0d0m+7mZM7nf6GlMBw+i+Y2zbP5mDENfnsKkjVDT0YKU9Ew8qligtjUr9fjLI41Gg6IoZGVl6ZyTy7k0Iue5sqZWrVocOHBAp2zKlClMmTJFu52VlcWuXbu0PwNcvXpVZxv+G67NKXN3d89V9rCsrCwURUGj0WBoaKjzXGF/l+gtUcbGxpKZmamdtZTDxcWF6OjoPPdp1aoVq1evJjAwkNTUVDIyMnjhhRf47LPP8n2d2bNnM3369Fzl27dvx6IIJ6RDQ0MLXVeULGmLsuNRbXEpy5zfghbw1U8fUOvuDdb/8C7vdhnJK18oKKhQoRBYIwtfl5K9Dq4iMDIywtXVlcTERNLT03M9n5CQoIeoyr709HRSUlLYu3cvGRkZOs8lJycX6hgqpaSv1MzHP//8Q7Vq1Th48CC+vr7a8g8//JDvvvuO8+fP59rn7NmzPPfcc4wZM4bOnTsTFRXFhAkT8PHxYcWKFXm+Tl49Sjc3N2JjY7GxsXlknBqNhtDQUDp16lRm1yGsLKQtyo7CtkVUXCrt5+/FMjWJhVvm0fFKGACf+/Zmfps3UFQGGKhg97i20rN8hNTUVK5fv46npydmZv99VoqikJCQgLW1tXYWqPhPamoqERERuLm56XxukJ0PHB0diYuLKzAf6K1H6ejoiKGhYa7eY0xMTK5eZo7Zs2fTunVrJkyYAGRfY2NpaUmbNm2YOXMmarU61z6mpqaYmprmKjc2Ni7SL9ui1hclR9qi7HhUW7g7GjO7RyMmbzzNkB7vMX7f9ww/vJ63D63D/X40E7oGk2Zkws24dNwdrUsx8vInMzMTlUqFgYGBzmUgOUOOOc8JXQYGBqhUqjy/q4X9PaK3T9XExARvb+9cQzehoaG0atUqz32Sk5NzfRFyxpz11DEWQjxCziSf1UNb03rdF0zoGozGwJAXzu1l9Y9TcEqJx9OxdK7LE+Jx6PXPj7Fjx7J8+XJWrlzJuXPnGDNmDJGRkQwbNgyASZMmERQUpK3fvXt3Nm7cyNKlS7l69SoHDhxg1KhRtGjRgqpVq+rrbQghHiFnVZ8mbvY0f38MA3p/QLypJc1vnmPXpsmo70QBEBWXIpeRPIJ0CoqmOD4vva7MExgYyJ07d5gxYwZRUVE0bNiQkJAQPDw8AIiKitK5prJ///4kJCTw+eefM27cOOzs7OjQoQMfffSRvt6CEKKIAn3cabssmKv929PwzdewjoyA1q3ZNu9r3jqdSZaCrBWbh5zRs/T0dMzNZbWjwsqZsPMkp2v0NplHX+Lj47G1tX3kydscGo2GkJAQunbtKufF9EzaouwotraIioIuXeDvv0kwtWBQz6kccWsIgKFKxf6JfrIE3r8URSEyMhKNRkPVqlW1p6GysrJITEzEyspKzlE+QFEUkpOTiYmJwc7OLs85LIXNB7LWqxBCf9Rq2LOHOP8AbMMO893a/zHixYnsqN2STEWRtWIfoFKpUKvVhIeHc+3aNW25oiikpKRgbm4us17zYGdn98R3FpFEKYTQLzs7kjf/xp9tn8f/0mGW/jyL0d0nsK1+G5nk8xATExNq166tcx2lRqNh7969tG3bVkZaHmJsbJxrkYHHIYlSCKF3alcH9n2zml8GD+LFs7v5bPNcwpq5oLbtqu/QyhwDAwOd6wENDQ3JyMjAzMxMEmUJkQFtIUSZ0Nu3Bi32beHWK69jqGTxzPSx3P90scyCFXoniVIIUWaoHaxw+fFbePttUBTsRr/N5hHTZTF1oVcy9CqEKFsMDIj64CN++/MGg8N+Zs62z8k0MGQyKtrWcZLJPaLUSaIUQpQ54XeSmek3CIOsLAYe28xHv39KhoEhEbEtJVGKUidDr0KIMsfL0RIDAxUzOg7h22bPY4DCvJCF1Nvzm75DE5WQJEohRJmjtjVndo9GGBoYMK3TUNY07YKhkoX90EHcXbdRJviIUiWJUghRJuUspv7Dm61ov2sDvP46ZGRg8fprfDJ1hUzwEaVGEqUQoszKWUxdbW9J1MIl7Kzpg1lGOis2zKB+9BUmbzwtPUtR4iRRCiHKhfC4dIa/OJE/3Rpik57MN+um4nbnBhGxhbtLvRCPSxKlEKJc8HK0RGNiyqCeUznlUhPH5Di+Wf8+NZREfYcmKjhJlEKIciFngk+KmSUDXnmfSFsXPO5H4dLnFUhK0nd4ogKT6yiFEOVGoI87bes4ERGbjFn/xtC5A4SFkdrzFY4vWomXq61cZymKnfQohRDlSs4EH2fvxrB5Mxmmppht+52rrw6k9eydMhNWFDtJlEKIcivqqWa83XUcWah44+Tv9Du6WWbCimIniVIIUW6FxyaxtU4rZrcfAMB7u1bQ5nKYzIQVxUoSpRCi3PJytMRABV+1eJm1jTphqGTx2eaPqBUTru/QRAUiiVIIUW49uNTde52Hc9itEdbpKTj1eQVu39Z3eKKCkFmvQohy7cGZsJ7DQ6BTe7hyhbQevTj+1Y94qu1kJqx4ItKjFEKUezkzYV29qsPmzWgsLDHdv5dzbwyTNWHFE5NEKYSoUKKqeTGy82gABh7bzIundslMWPFEJFEKISqUnJmwn/oGAjB72+fUj7okM2HFY5NEKYSoUHJmwn7S5nXt3Ua+2PQhXgbSoxSPRxKlEKJCyZkJa2BgyJhu44iwr0q1+Nu4vj2UqHtJctNnUWQy61UIUeE8OBPWsudP4O8HISGsfmk4n/v2xkAFs3s0ItDHXd+hinJAepRCiAopZyasU5tnuD93AQBj9n2P77W/yVKQCT6i0CRRCiEqvLMBvVjf8DkMlSw+3TIXp8S7ZCqKTPARhSKJUghR4Xk5WTGt8zDOOXnilHSfT36djxEKno4W+g5NlAOSKIUQFZ7a1pxpgT6MfGkSycamPHvtLzYm7JMVe0ShSKIUQlQKgT7ufPdxX/754GMAGn8xHw4d0nNUojyQRCmEqDTUtubUeudteO01yMzM/vf+fX2HJco4SZRCiMpFpYJly8DLC65dI2XAIA5evi0zYEW+JFEKISofGxv48UeyDI0w/3kjG4Jny+LpIl96T5RLlizBy8sLMzMzvL292bdvX4H109LSmDJlCh4eHpiamlKzZk1WrlxZStEKISqKqLqNWNC6DwDTQ5dR9f4tubZS5EmviXLt2rUEBwczZcoUTpw4QZs2bQgICCAyMv+/6nr37s3OnTtZsWIFFy5cYM2aNdSrV68UoxZCVAThsUksadmTsGoNsE5PYf6vC1AyM+TaSpGLXhPlggULGDRoEIMHD6Z+/fosXLgQNzc3li5dmmf9rVu3smfPHkJCQnjuuefw9PSkRYsWtGrVqpQjF0KUd16OlmBoyJhuY0k0MafljTMMO7JRrq0UuegtUaanp3Ps2DH8/f11yv39/Tl48GCe+2zevJnmzZszd+5cqlWrRp06dRg/fjwpKTJUIoQompzF06Ps1Uzv+CYA4w78gPrqeT1HJsoavS2KHhsbS2ZmJi4uLjrlLi4uREdH57nP1atX2b9/P2ZmZmzatInY2FiGDx/O3bt38z1PmZaWRlpamnY7Pj4eAI1Gg0ajeWScOXUKU1eULGmLsqOitEWPpmp8veyJvPM0yQZXsfhtC0pQEBmHDoGpqb7DK5SK0hb6UNjPTO93D1GpVDrbiqLkKsuRlZWFSqVi9erV2NraAtnDt7169WLx4sWYm+deZWP27NlMnz49V/n27duxsCj8EEtoaGih64qSJW1RdlSkttjTuxcd9u/F9PRprg4cyPnXX9d3SEVSkdqitCQnF+58tN4SpaOjI4aGhrl6jzExMbl6mTnUajXVqlXTJkmA+vXroygKN27coHbt2rn2mTRpEmPHjtVux8fH4+bmhr+/PzY2No+MU6PREBoaSqdOnTA2Ni7s2xMlQNqi7KiobaEyNYVXX6XOxo0k9x5IlfatUdua6TusAlXUtigNOSOMj6K3RGliYoK3tzehoaG8/PLL2vLQ0FBefPHFPPdp3bo169evJzExESsrKwAuXryIgYEB1atXz3MfU1NTTPMYQjE2Ni7Sl6qo9UXJkbYoOypcWwQGcm3lD3hs34zVsDfxH7CI6b29y8V9KytcW5SCwn5eep31OnbsWJYvX87KlSs5d+4cY8aMITIykmHDhgHZvcGgoCBt/T59+lClShUGDBjA2bNn2bt3LxMmTGDgwIF5DrsKIURRRMWl8HL9V7ltaUftO9cJ3rdarq0U+k2UgYGBLFy4kBkzZtC0aVP27t1LSEgIHh4eAERFRelcU2llZUVoaCj379+nefPmvP7663Tv3p1PP/1UX29BCFGBhMcmcdfMhsmd3wZgyJFNNLx5Xq6trOT0Ppln+PDhDB8+PM/nVq1alausXr16ctJaCFEivBwtMVBBaO1n2NSgPS+f3c3crZ9iM3+QvkMTeqT3JeyEEKKsyLm20lClYkbHIdwxt6Hu7Wuoly7Sd2hCj/TeoxRCiLIk0MedtnWciIhNxqDpZzCoH8ycCT17wlNP6Ts8oQfSoxRCiIeobc3xrVkF+wF9oXt30GhIHzCQgxdvycSeSkgSpRBC5EelgiVL0FhaYRJ2hNDh/5PbcVVCjzX0qigKO3bs4ODBg0RHR6NSqXBxcaF169Z07Ngx35V1hBCivImyrsLnz/bjw22LGb/3O7bWacXkjadpW8cJta1cllYZFLlHefPmTZ5++mkCAgLYtGkTV69e5fLly2zatIkuXbrQvHlzbt68WRKxCiFEqQuPTeKHJp0Jq9YAS00q03d8QaaiyCUjlUiRE+Xw4cNxcHDg+vXrnDx5km3btrF9+3ZOnjzJ9evXsbOzY8SIESURqxBClDovR0tUBgZM7jwCjYEh/pcO0+XSYbkdVyVS5ES5c+dOFixYgFqtzvWcWq1m3rx57Nixo1iCE0IIfcu5ZOSqsydftOwJwIL9y1EbZOg5MlFaipwozc3NuXv3br7P37t3T5aTE0JUKIE+7uyf6IfPV/PJ8PTCIiYa/vc/fYclSkmRE+Wrr75Kv3792LBhA3FxcdryuLg4NmzYwIABA+jTp0+xBimEEPqmtjWn5VPVMfpiWXbBZ5/BsWP6DUqUiiLPep0/fz4ZGRm8/vrrZGRkYGJiAkB6ejpGRkYMGjSIjz/+uNgDFUKIMsHfH157DdasIX3oMI6u+Q0vZ2uZAVuBFTlRmpiYsHTpUj766COOHj3KrVu3AHB1dcXb27tQ93gUQohybf58NJu3YHLsKFtGfsDaZl2Y3aNRubgdlyi6x17CzsbGhg4dOhRnLEIIUS5EWdjx1TOvMXXnV7y7ZxXb6vjKtZUVWJETZWFvaTVq1KgiByOEEOVBeGwS3zzdjV6ndtAgJpx3d6/i3a6jiYhNlkRZARU5UX7yySc629evX0etVmNk9N+hVCqVJEohRIXl5WiJYmjIe52Gs3H1BAJPhbKhiT+ejjLKVhEVOVGGh4frbFtbW7Nnzx5q1KhRbEEJIURZlnNt5eSNKtY1eo7ep3bw5ZGvsbcao+/QRAmQRdGFEOIx5Fxb6bX8c7Js7bC/eBa+/FLfYYkSIIlSCCEek9rWHJ8WdTH4YEZ2wXvvwZ07+g1KFDtJlEII8aTeegsaNYK7d2XFngqoyIkyPj5e56FSqUhMTMxVLoQQlYaRUfZKPYDyxRec3PKH3OC5AilyorSzs8Pe3l77SExMpFmzZtrtnOeFEKJSadeOyE7dUWVlkT78bVrP3ik3eK4gijzr9Y8//iiJOIQQolyLikvhtdo9CN29nRY3ztLt7B4mqwxkEYIKoMiJMi0tDT8/P4yNjUsiHiGEKJfCY5O4ae3E4md6M2Hfd0zcvYrQWs/IIgQVQJGHXocNG4aTkxOBgYH88MMP3L9/vwTCEkKI8sXL0RIDFSz3eYkbNs5UTYjlrSM/yQ2eK4AiJ8qrV6+yd+9eGjVqxMKFC3F1daVjx458+umnRERElECIQghR9uUsQpBhYsaHfgMBGH50I+q423qOTDypx7o8pHHjxrz33nscOXKEK1eu8Morr7B161bq169PkyZNmDp1KkePHi3uWIUQokzLWYQgaN440lo9i1FaGrz7rr7DEk/oia+jrFatGsOGDSMkJITY2FimTp1KREQEXbp0YdasWcURoxBClBtqW3N8azliuvgzUKngxx9h/359hyWewGPfZgtg586d7Ny5k5iYGLKysv47qJERMTEx3JEVKoQQlVXTpjB4MHz1FQQHw5EjYCBrvJRHj91q06dPx9/fn507dxIbG8u9e/e0j/v372NgYICTk1NxxiqEEOXLzJlgYwPHjnFpwVJZhKCceuwe5bJly1i1ahV9+/YtzniEEKLicHbmr6ARNPl8NtbTp9IxypWpr/oQ6OOu78hEETx2jzI9PZ1WrVoVZyxCCFGhRMWlEGjmQ6StC66Jdxny509M3nhaepblzGMnysGDB/PDDz8UZyxCCFGhhMcmkWpowpz2AwAY+udGHONvExGbrOfIRFEUaeh17Nix2p+zsrL48ssv2bFjB40bN861Us+CBQuKJ0IhhCinchYhCKnbmiPVG9Dixlne3fsdnrN66zs0UQRFSpQnTpzQ2W7atCkAp0+f1ilXqVRPFpUQQlQAOYsQTN54mpkdBrP527H0OL0TLp2B5s31HZ4opCIlSlkQXQghiibQx522dZyIiG1JcmoYFuvWwLhxsHt39nWWoszT+0U9S5YswcvLCzMzM7y9vdm3b1+h9jtw4ABGRkbaXq0QQpRValtzfGtWwWLeR2BmBnv3cnfNeg5eiZWJPeWAXhPl2rVrCQ4OZsqUKZw4cYI2bdoQEBBAZGTB93CLi4sjKCiIjh07llKkQghRDNzc4N+5HvdHjiVo2QFaz9kl960s4/SaKBcsWMCgQYMYPHgw9evXZ+HChbi5ubF06dIC9xs6dCh9+vTB19e3lCIVQojiET1sFLEWdtS4e5M+J38nS0EuGSnj9JYo09PTOXbsGP7+/jrl/v7+HDx4MN/9vv76a65cucK0adNKOkQhhCh2V9MN+eTZPgAEH1iDTWoimYoil4yUYU+01uuTiI2NJTMzExcXF51yFxcXoqOj89zn0qVLTJw4kX379mFkVLjQ09LSSEtL027Hx8cDoNFo0Gg0j9w/p05h6oqSJW1RdkhbPL7qtqasa9qZAUc3U+vuDYYfXs9cvwFUszV5rM9T2uLxFfYz01uizPHwpSSKouR5eUlmZiZ9+vRh+vTp1KlTp9DHnz17NtOnT89Vvn37diwsCn9D1dDQ0ELXFSVL2qLskLZ4PD1rqpjtN4AVP33AgKO/EPdCZ04c2MWJR++aL2mLoktOLlwvXqUoilLCseQpPT0dCwsL1q9fz8svv6wtHz16NCdPnmTPnj069e/fv4+9vT2GhobasqysLBRFwdDQkO3bt9OhQ4dcr5NXj9LNzY3Y2FhsbGweGadGoyE0NJROnTrlWlRBlC5pi7JD2uLJRd1PwaZ7V+z+PEDWq6+S+e23j3UcaYvHFx8fj6OjI3FxcQXmA731KE1MTPD29iY0NFQnUYaGhvLiiy/mqm9jY8OpU6d0ypYsWcKuXbvYsGEDXl5eeb6OqakppqamucqNjY2L9KUqan1RcqQtyg5pi8fn7mQMixdB8+YY/PgjBuPGPdEiBNIWRVfYz0uvQ69jx46lb9++NG/eHF9fX7788ksiIyMZNmwYAJMmTeLmzZt8++23GBgY0LBhQ539nZ2dMTMzy1UuhBDlgrc3vPEGfP89TJgAu3bJIgRlkF4TZWBgIHfu3GHGjBlERUXRsGFDQkJC8PDwACAqKuqR11QKIUS5NnMmrF+fvVLPb79Bt276jkg8RO8r8wwfPpyIiAjS0tI4duwYbdu21T63atUqdu/ene++77//PidPniz5IIUQoqR4eMDo0QBoxk/g4IVouaayjNF7ohRCiEpv0iTSbO0xvnCeLaM/lNV6yhhJlEIIoWdRKlPm+LwCwNj932OWliKr9ZQhkiiFEELPwmOT+L5pABF2apyS7jPkyCZZracMkUQphBB65uVoSaaRMXPb9QPgzSMbcUm6h6dj4RdFESVHEqUQQuhZzg2et9V7lhPqulhqUvnh5u+obc31HZpAEqUQQpQJgT7u7J/UAeNP5gNQ85cf4dw5PUclQBKlEEKUGWpbcxoGPg8vvQRZWfDuu/oOSSCJUgghyp45c8DQELZs4c5v2zl4JVZmwOqRJEohhChr6taFN98E4Obgt3n9y0NybaUeSaIUQogy6NaYd0g0Madx9CW6n9tHloJcW6knkiiFEKIMumJgxbKWPQGYsPdbTDI0cm2lnkiiFEKIMsjL0ZKVLV4i2soBt7hbBB3fgqFKJddW6oEkSiGEKIPUtuZMC/ThkzZ9ARh5cC3zOlaXayv1QBKlEEKUUYE+7gR/P5Okug2wTUvi5d+/0XdIlZIkSiGEKMPUDlZYLsxehIDPPoOrV/UbUCUkiVIIIcq6zp2hUyfQaGDSJH1HU+lIohRCiLJOpYJ587L/XbcODh3Sd0SViiRKIYQoDxo3hgEDAEgPHsPBy7flmspSIolSCCHKiw8+IMPMHJMjf/Ld+AW0nrOL9cdu6DuqCk8SpRBClBNRlvYs9n4JgIl7VmGk0fDeL2e5n6bfuCo6SZRCCFFOhMcm8UWLHtyycsDjfjR9j/9KlgK3U1X6Dq1Ck0QphBDlhJejJamm5sxr8wYAow/+iGNKHE5mip4jq9gkUQohRDmhtjVndo9G/NzoOc4418AmLYkfbvyOnam+I6vYJFEKIUQ5Eujjzt7Jz8GCBQDU/mUN1pFy+62SJIlSCCHKGbWtOU+9/iK8/DKqzEye+vprfYdUoUmiFEKI8urjj1GMjXE5cQLV1q36jqbCkkQphBDlVc2aZI0cCYAybjyHzkfJIgQlQBKlEEKUY1mTJpFobYvxpYv8PnI6refsYm2YnLMsTpIohRCiHIvClFmtXgdgzP4fsE5OYPLG09KzLEaSKIUQohy7dieZtY07c87JE/vUBIIP/ECmohARm6zv0CoMSZRCCFGOeVSxIMvAgA86DAag7/HfqHPnOp6OFnqOrOKQRCmEEOWY2taMwBpZHPZqSmitlhgpWaw6tw61rbm+Q6swJFEKIUQ55+uisHtcWxyXfYpibEzVA3/A77/rO6wKQxKlEEJUAGpbM5p1bIFq1KjsgpEjITVVv0FVEJIohRCiIpk6FdRquHIFPv5Y39FUCJIohRCiIrGx0a4Dq8yaxbHdx+RSkSek90S5ZMkSvLy8MDMzw9vbm3379uVbd+PGjXTq1AknJydsbGzw9fVl27ZtpRitEEKUA4GB3GreClVqKvcGvyWLEDwhvSbKtWvXEhwczJQpUzhx4gRt2rQhICCAyHxWwt+7dy+dOnUiJCSEY8eO4efnR/fu3Tlx4kQpRy6EEGVXVHwqrzfrS7qBEc9dCaPDxT9lEYInoNdEuWDBAgYNGsTgwYOpX78+CxcuxM3NjaVLl+ZZf+HChbzzzjv4+PhQu3ZtZs2aRe3atdmyZUspRy6EEGVXeGwSlx3cWN7iJQDe3/EFJmkpsgjBYzLS1wunp6dz7NgxJk6cqFPu7+/PwYMHC3WMrKwsEhIScHBwyLdOWloaaWlp2u34+HgANBoNGo3mka+RU6cwdUXJkrYoO6Qtyo682qK6rSkGKvjM91VeOLuH6vExjDnwA9Um+0ubPaCwn4XeEmVsbCyZmZm4uLjolLu4uBAdHV2oY8yfP5+kpCR69+6db53Zs2czffr0XOXbt2/HwqLwK1eEhoYWuq4oWdIWZYe0RdnxcFv09lKx9qop//MfztcbpjPo6M/s/W45J2rU0FOEZU9ycuF62HpLlDlUKpXOtqIoucrysmbNGt5//31++eUXnJ2d8603adIkxo4dq92Oj4/Hzc0Nf39/bGxsHvk6Go2G0NBQOnXqhLGx8SPri5IjbVF2SFuUHfm1RVdgeFwqkXd9SEr/G8vNm2j3ww9k7tsHhob6C7gMyRlhfBS9JUpHR0cMDQ1z9R5jYmJy9TIftnbtWgYNGsT69et57rnnCqxramqKqalprnJjY+Mi/Qcvan1RcqQtyg5pi7Ijr7ZwdzTG3dEaln4Ou3dicPQoCZ8u4UzPILwcLSv9MneF/e7qbTKPiYkJ3t7euYYLQkNDadWqVb77rVmzhv79+/PDDz/w/PPPl3SYQghR/lWtCnPmAGD4v/cYO/83uWSkCPQ663Xs2LEsX76clStXcu7cOcaMGUNkZCTDhg0DsodNg4KCtPXXrFlDUFAQ8+fP55lnniE6Opro6Gji4uL09RaEEKJciHo1iGPV6mGVnsLM7YvJylLkkpFC0muiDAwMZOHChcyYMYOmTZuyd+9eQkJC8PDwACAqKkrnmsovvviCjIwMRowYgVqt1j5Gjx6tr7cghBDlQvjdFCZ2HkmaYfa1lS+d3S33rSwkvU/mGT58OMOHD8/zuVWrVuls7969u+QDEkKICsjL0ZIrzh582uo1Juz7jvd3fMEhz6Zy38pC0PsSdkIIIUqe2tac2T0asfyZXpx2qYldaiI/nf4etY2ZvkMr8yRRCiFEJRHo487uKZ3IXL4CxciI6ru3wbp1+g6rzJNEKYQQlYja1pwm3dqheu+97IK33+bWxXAOXomViT35kEQphBCV0aRJ0LQpxMZytmtv+nx5WC4ZyYckSiGEqIxMTLi9dDlphsb4XTnKGydCyFKQS0byIIlSCCEqqUtOHsxuPwCA9/5YQc3Y63LJSB4kUQohRCXl5WjJd827sdezGWYZ6Sz8dR5mWRlyychDJFEKIUQlpbY1Z1bPJrz7/BjumtvQ6NYVNkX9XunXgH2YJEohhKjEAn3c2fhhL27P/wyA+t9/Ab/+SlRcisyE/ZfeV+YRQgihX2pbc9RvBcGF47BoEWmvv0HvPgu4buuCgQpm92hEoI+7vsPUG+lRCiGEyDZ3LunezTGNj+Pzn+dgkqGRmbBIohRCCJHDxIRTC77kvpkVTaIvMfmPFQCVfiasJEohhBBaVZvUY1y3cQD0P/4rvU7twFClqtQzYSVRCiGE0FLbmuP/zkA+a/UqAB9u+5wva6ZU6pmwkiiFEELoCPRxp9evy7nTuTummRl0nPwWREToOyy9kUQphBAiF7W9JVV+WgPNmsHt2/DCC5CQoO+w9EISpRBCiLxZWsIvv4CLC5w6RWqvVzh0PqrSzYCVRCmEECJ/bm7wyy9kmJphtn0bt3q+xrOzd1Squ4xIohRCCFGgqHqNefOFd9EYGPLS2T38L/RLJv90qtL0LCVRCiGEKFB4bBK7avgw7vkxZKGi//FfGbl/NRGxyZViqTtZwk4IIUSBvBwtMVDB5gbtsU1N5IPQZQQfWMOROWpaO3YgS6FCL3UnPUohhBAFUtuaM7tHIwxVKr57uhvz2gYB0GL5Asbt/gYUpUIvdSc9SiGEEI8U6ONO2zpORMQm4zmpAxFzauI5ZxojDq/HQpPK9I5vkglExCZXuMUJpEcphBCiUNS25vjWrILa1hzTiRP4n/9wAAYc28Lc3xdhlpWBhYlBhTtnKT1KIYQQRaa2NafhzHeZYGzKnJBF9D61g6cz7xGYNpY75rYV6pyl9CiFEEI8lkAfd8au/pCLX60m08qKWmePsfHbcdSKjaxQ5ywlUQohhHhsaltz6g96lb/X/U6krQse96PZ+N14As7vrzC355JEKYQQ4om5tvLm5X4L+LP6U9ikJ7P0lznM/X0RXmZZ+g7tiUmiFEII8cTUtua888azBL02i8XPvEIWKnr/HYqrX2ti/9hXrif4SKIUQghRLAJ93Nk9pRNPf7+Ue1t+h+rV4dIl7Du25/Ir/Ql4f3O5XCNWEqUQQohik3MJSZVunYned5hfGrTDUMki6MRv7PziTY5PW0DUvaRytfSdJEohhBAl4mqmKaO7T+C1V2dxsYo7VVLi+ej3RRg1bcrsoGm88cVBWs/Zpe1lltXkKddRCiGEKBE5a8Qe8mhM1wGf0u/YFkYf/BGnyMt8Gvkxox3WsOSZV5iuSed+ioaPfj9fJteNlR6lEEKIEvHgGrEZhkasatmDr77dybw2b3DfzIqad28wP+QTDn3eF9OxY6gdEwGgcw1mXr3M0u55So9SCCFEidFZI9bRAoDWJ15llfcL9D3xG6+f+J3q8TH0P7aF/se2cMHRne21fQmt3ZKv93qw/GCETi8TYNLGU6Xa89R7j3LJkiV4eXlhZmaGt7c3+/btK7D+nj178Pb2xszMjBo1arBs2bJSilQIIcTjeHCN2JxeZoqZJUufeQW/Ycv5be5Kfq/bCo2BIXVjIxl5aC2bvx3L4Nfa8OmmOfQ7toV60Vd5b/1JbZIESm31H732KNeuXUtwcDBLliyhdevWfPHFFwQEBHD27Fnc3XP/hRAeHk7Xrl0ZMmQI33//PQcOHGD48OE4OTnRs2dPPbwDIYQQRfVwL1Nta87a9h1pufoQba+E4X/pMB0jT+CcdI9uF/bT7cJ+ANIMjbnk6M4FJw/OO3qys1YLrlapXuJ3LNFrolywYAGDBg1i8ODBACxcuJBt27axdOlSZs+enav+smXLcHd3Z+HChQDUr1+fo0ePMm/ePEmUQghRjuT0LnP8lzyfw9PRgrupqQSP/ZLm18/Q4sYZnr55Duv0FBreukLDW1cAuG1lzzVHN+2QbknRW6JMT0/n2LFjTJw4Uafc39+fgwcP5rnPoUOH8Pf31ynr3LkzK1asQKPRYGxsnGuftLQ00tLStNtxcXEA3L17F41G88g4NRoNycnJ3LlzJ8/ji9IjbVF2SFuUHRWpLUyAOnZARjIYQaeBAcz8rQafeXfHkCxmPW2N3dWLnNh2iFq3IzlZpRqTnquOSUYyd+4UfU3ZhIQEABRFKbCe3hJlbGwsmZmZuLi46JS7uLgQHR2d5z7R0dF51s/IyCA2Nha1Wp1rn9mzZzN9+vRc5V5eXk8QvRBCiNLW5+GCy4fZ8w0MecLjJiQkYGtrm+/zep/1qlKpdLYVRclV9qj6eZXnmDRpEmPHjtVuZ2VlcffuXapUqUKLFi0ICwvTqe/j46NTFh8fj5ubG9evX8fGxqZwb6qYPBxLaRyjsPUfVa+g5/N6rjBlla0tCrtPcbdFfuUPlpX3tnic4xTH/w1pi+I5TnH9nmrevDm7du2iatWqBR5Hb4nS0dERQ0PDXL3HmJiYXL3GHK6urnnWNzIyokqVKnnuY2pqiqmpqU6ZnZ0dAIaGhrm+WHmVAdjY2JT6lzC/WEryGIWt/6h6BT1f2M+9srdFYfcp7rbIrzyvsvLaFo9znOL4vyFtUTzHKa7fU0ZGRlSvXv2Rx9Hb5SEmJiZ4e3sTGhqqUx4aGkqrVq3y3MfX1zdX/e3bt9O8efPHGpsfMWJEocr0pThiKeoxClv/UfUKer6wn3tlb4vC7lPcbZFfeVlpj+KKQx//N6Qtiuc4pfF7SoeiRz/++KNibGysrFixQjl79qwSHBysWFpaKhEREYqiKMrEiROVvn37autfvXpVsbCwUMaMGaOcPXtWWbFihWJsbKxs2LChxGKMi4tTACUuLq7EXkMUjrRF2SFtUXZIW5Q8vZ6jDAwM5M6dO8yYMYOoqCgaNmxISEgIHh4eAERFRREZ+d8tWby8vAgJCWHMmDEsXryYqlWr8umnn5bopSGmpqZMmzYt1/CtKH3SFmWHtEXZIW1R8lSK8oh5sUIIIUQlpvcl7IQQQoiyTBKlEEIIUQBJlEIIIUQBJFEKIYQQBZBEKYQQQhRAEmUxCg8Px8/PjwYNGtCoUSOSkpL0HVKlZWRkRNOmTWnatKn27jRCf5KTk/Hw8GD8+PH6DqVSS0hIwMfHh6ZNm9KoUSO++uorfYdULsjlIcWoXbt2zJw5kzZt2nD37l1sbGwwMtL7crqVkqOjI7GxsfoOQ/xrypQpXLp0CXd3d+bNm6fvcCqtzMxM0tLSsLCwIDk5mYYNGxIWFpbvEqAim/Qoi8mZM2cwNjamTZs2ADg4OEiSFAK4dOkS58+fp2vXrvoOpdIzNDTEwiL73o2pqalkZmY+8hZTohIlyr1799K9e3eqVq2KSqXi559/zlVnyZIleHl5YWZmhre3N/v27Sv08S9duoSVlRUvvPACTz/9NLNmzSrG6CuWkm4LyL6jgre3N88++yx79uwppsgrntJoi/Hjx+d5I3aRW2m0x/3792nSpAnVq1fnnXfewdHRsZiir7gqTZcnKSmJJk2aMGDAgDyXvFu7di3BwcEsWbKE1q1b88UXXxAQEMDZs2dxd3cHwNvbW+cm0Dm2b9+ORqNh3759nDx5EmdnZ7p06YKPjw+dOnUq8fdW3pR0W1StWpWIiAiqVq3K6dOnef755zl16lSp31mhPCjptggLC6NOnTrUqVMn3xuyi/+Uxv8NOzs7/vrrL27dukWPHj3o1atXvndsEv/S60qzegIomzZt0ilr0aKFMmzYMJ2yevXqKRMnTizUMQ8ePKh07txZuz137lxl7ty5TxxrRVcSbfGwLl26KGFhYY8bYqVREm0xceJEpXr16oqHh4dSpUoVxcbGRpk+fXpxhVyhlcb/jWHDhinr1q173BArjUoz9FqQ9PR0jh07hr+/v065v79/of8K9vHx4datW9y7d4+srCz27t1L/fr1SyLcCq042uLevXvav6hv3LjB2bNnqVGjRrHHWtEVR1vMnj2b69evExERwbx58xgyZAhTp04tiXArvOJoj1u3bhEfHw9kn57Yu3cvdevWLfZYK5pKM/RakNjYWDIzM3MNP7i4uOS6UXR+jIyMmDVrFm3btkVRFPz9/enWrVtJhFuhFUdbnDt3jqFDh2JgYIBKpWLRokU4ODiURLgVWnG0hSg+xdEeN27cYNCgQSiKgqIovP322zRu3Lgkwq1QJFE+QKVS6WwripKrrCABAQEEBAQUd1iV0pO0RatWrTh16lRJhFUpPen/ixz9+/cvpogqtydpD29vb06ePFkCUVVsMvRK9jV3hoaGuf4qi4mJkZPcpUzaouyQtihbpD30RxIlYGJigre3N6GhoTrloaGhtGrVSk9RVU7SFmWHtEXZIu2hP5Vm6DUxMZHLly9rt8PDwzl58iQODg64u7szduxY+vbtS/PmzfH19eXLL78kMjKSYcOG6THqiknaouyQtihbpD3KKD3OuC1Vf/zxhwLkevTr109bZ/HixYqHh4diYmKiPP3008qePXv0F3AFJm1RdkhblC3SHmWTrPUqhBBCFEDOUQohhBAFkEQphBBCFEASpRBCCFEASZRCCCFEASRRCiGEEAWQRCmEEEIUQBKlEEIIUQBJlEIIIUQBJFEKUcbt3r0blUrF/fv3S/21VSoVKpUKOzu7Auu9//77NG3aVGc7Z9+FCxeWaIxClDRJlEKUIe3btyc4OFinrFWrVkRFRWFra6uXmL7++msuXrxYpH3Gjx9PVFQU1atXL6GohCg9lWZRdCHKKxMTE1xdXfX2+nZ2djg7OxdpHysrK6ysrDA0NCyhqIQoPdKjFKKM6N+/P3v27GHRokXaYcuIiIhcQ6+rVq3Czs6OX3/9lbp162JhYUGvXr1ISkrim2++wdPTE3t7e0aOHElmZqb2+Onp6bzzzjtUq1YNS0tLWrZsye7dux8r1jlz5uDi4oK1tTWDBg0iNTW1GD4BIcomSZRClBGLFi3C19eXIUOGEBUVRVRUFG5ubnnWTU5O5tNPP+XHH39k69at7N69mx49ehASEkJISAjfffcdX375JRs2bNDuM2DAAA4cOMCPP/7I33//zSuvvEKXLl24dOlSkeJct24d06ZN48MPP+To0aOo1WqWLFnyRO9diLJMhl6FKCNsbW0xMTHBwsLikUOtGo2GpUuXUrNmTQB69erFd999x61bt7CysqJBgwb4+fnxxx9/EBgYyJUrV1izZg03btygatWqQPZ5xK1bt/L1118za9asQse5cOFCBg4cyODBgwGYOXMmO3bskF6lqLCkRylEOWRhYaFNkgAuLi54enpiZWWlUxYTEwPA8ePHURSFOnXqaM8fWllZsWfPHq5cuVKk1z537hy+vr46ZQ9vC1GRSI9SiHLI2NhYZ1ulUuVZlpWVBUBWVhaGhoYcO3Ys1wSbB5OrECI3SZRClCEmJiY6E3CKS7NmzcjMzCQmJoY2bdo80bHq16/P4cOHCQoK0pYdPnz4SUMUosySRClEGeLp6cmff/5JREQEVlZWODg4FMtx69Spw+uvv05QUBDz58+nWbNmxMbGsmvXLho1akTXrl0LfazRo0fTr18/mjdvzrPPPsvq1as5c+YMNWrUKJZYhShr5BylEGXI+PHjMTQ0pEGDBjg5OREZGVlsx/76668JCgpi3Lhx1K1blxdeeIE///wz35m1+QkMDGTq1Km8++67eHt7c+3aNd56661ii1OIskalKIqi7yCEEGWTSqVi06ZNvPTSS4+1v6enJ8HBwblWGxKiPJEepRCiQK+99lqRl6KbNWsWVlZWxdojFkJfpEcphMjX5cuXATA0NMTLy6vQ+929e5e7d+8C4OTkpLd1aoUoDpIohRBCiALI0KsQQghRAEmUQgghRAEkUQohhBAFkEQphBBCFEASpRBCCFEASZRCCCFEASRRCiGEEAWQRCmEEEIUQBKlEEIIUYD/A9fTUEFXe4kEAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "tm = np.logspace(np.log10(to[0]), np.log10(to[-1]), 100)\n", "hm = ml.head(0, 0, tm, layers=1)\n", "plt.semilogx(to, ho / H0, \".\", label=\"obs\")\n", "plt.semilogx(tm, hm[-1] / H0, \"r\", label=\"ttim\")\n", "plt.ylim([0, 1])\n", "plt.xlabel(\"time [d]\")\n", "plt.ylabel(\"h/H0\")\n", "plt.title(\"Model Results - Three layers\")\n", "plt.legend()\n", "plt.grid()" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[6.049085184433308, 0.00021324705409226575, 0.002855742798130468]" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ca.parameters[\"optimal\"].values.tolist() + [ca.rmse()]" ] }, { "cell_type": "code", "execution_count": 11, "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", "
Comparison of parameter values and error under different models
 k [m/d]Ss [1/m]RMSE[m]
AQTESOLV4.0340000.0003830.002976
ttim6.0490850.0002130.002856
\n" ], "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ta = pd.DataFrame(\n", " columns=[\"k [m/d]\", \"Ss [1/m]\", \"RMSE[m]\"],\n", " index=[\"AQTESOLV\", \"ttim\"],\n", ")\n", "ta.loc[\"ttim\"] = ca.parameters[\"optimal\"].values.tolist() + [ca.rmse()]\n", "ta.loc[\"AQTESOLV\"] = [4.034, 3.834e-04, 0.002976]\n", "ta.style.set_caption(\"Comparison of parameter values and error under different models\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Both models have similar (and very lowe) RMSE, but the hydraulic conductivity varies by 50%" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## References\n", "* Butler, J.J., Jr., 1998. The Design, Performance, and Analysis of Slug Tests, Lewis Publishers, Boca Raton, Florida, 252p.\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 }