{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 4. Slug Test - Dawsonville confined aquifer example\n", "\n", "**This test is taken from example of MLU.**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Import packages" ] }, { "cell_type": "code", "execution_count": 14, "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": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "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 MLU model (Carlson & Randall, 2012).\n", "\n", "This Slug Test was reported in Cooper Jr et al. (1967), and it was performed in Dawsonville, Georgia, USA. A fully penetrated well (Ln-2) is screened in a confined aquifer, located between depths 24 and 122 (98 m thick).\n", "\n", "The volume of the slug is 10.16 litres. Head change has been recorded at the slug well. Both the well and the casing radii of the slug well is 0.076 m.\n", "\n", "The conceptual model can be seen in the figure below:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "jupyter": { "source_hidden": true }, "tags": [ "hide-input" ] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAE6CAYAAAAlRjrfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABR7ElEQVR4nO3dd1gUV/s38O+CsLssRXpRimgUEBUBC5YgWLDxqPnFhgXUYAxiwxJbEIxIoth9LERFE2vsirFgi5qgoqISa2JQUEAQERSkn/cPX+ZxXdBd2GUo9+e69rqYM2fOuWdh92bOnJkRMMYYCCGEEMIbNb4DIIQQQuo7SsaEEEIIzygZE0IIITyjZEwIIYTwjJIxIYQQwjNKxoQQQgjPKBkTQgghPKNkTAghhPCMkjEhhBDCM0rGNdTt27cxZswYNGnSBCKRCNra2nB2dsaSJUvw8uVLvsNTupSUFISEhODmzZvV0l+3bt3QrVs3ueoJBALY2tqivJvVXbhwAQKBAAKBAFu3blVafFu3boVAIMDjx48V3jYkJAQCgUBpsQDA48ePuf0UCATQ0NCAoaEh2rVrh2nTpuHOnTtK7a+2K3u/3v+bKO936ufnBxsbG6X1W9ZHRa/z588rrS8+2djYwM/Pj+8wlKoB3wEQWT/99BMCAgLQokULzJw5Ew4ODigqKsK1a9ewYcMGxMbG4uDBg3yHqVQpKSkIDQ2FjY0NnJyc+A5Hio6ODhITE3H27Fl0795dat2WLVugq6uLnJwcnqKrXpMmTYKPjw9KS0vx6tUrxMfHY8uWLVizZg3Cw8Mxc+ZMvkOsEczNzREbG4umTZvy0n9UVBTs7Oxkyh0cHHiIhsiDknENExsbi2+++QY9e/bEoUOHIBQKuXU9e/bE9OnTceLECR4jrH+srKygo6ODLVu2SCXj169fY+/evRgxYgR++uknHiOsPlZWVujYsSO33LdvXwQFBeGLL77ArFmz4OjoiD59+vAYYc0gFAql3qfq5ujoCFdXV976J4qjYeoaZvHixRAIBIiMjJRKxGU0NTXxn//8h1suLS3FkiVLYGdnB6FQCBMTE4wePRpPnz6V2q5bt25wdHREXFwcunbtCi0tLdja2uKHH35AaWmpVN1Xr15h+vTpsLW15drs27cv7t+/z9UpLCzEokWLuH6NjY0xZswYZGRkSLVlY2OD/v374+DBg2jdujVEIhFsbW2xevVqrs758+fRrl07AMCYMWO4IbWQkBAu9vKGlMsb4gsNDUWHDh1gYGAAXV1dODs7Y/PmzeUOMSti7NixOHDgAF69esWV7d69GwAwbNiwcre5dOkSunfvDh0dHWhpaaFTp044duyYTL3Lly+jc+fOEIlEsLCwwJw5c1BUVFRum3v27IGbmxskEgm0tbXh5eWF+Pj4Ku1bVYnFYmzevBkaGhpYunQpV56RkYGAgAA4ODhAW1sbJiYm8PT0xMWLF6W2b9euHfr16ydV1qpVKwgEAsTFxXFlBw4cgEAgQEJCAtf++PHjYWlpyf0Ndu7cGadPn5Zqa8uWLWjTpg1EIhEMDAwwaNAg3Lt3T6qOn58ftLW18c8//6Bv377Q1taGpaUlpk+fjoKCAgBAUVERTExMMGrUKJn34NWrVxCLxQgKCgJQ/jC1vBhjWLduHZycnCAWi6Gvr48vv/wS//77r8JtVWT37t0QCARYu3atVPmCBQugrq6OmJgYrkzez1TZZz06Ohpt27aFWCyGvb09oqOjAbwbQre3t4dEIkH79u1x7do1qe3Lfgd37txB9+7dIZFIYGxsjMDAQOTl5X1yn3JycjBjxgw0adIEmpqaaNSoEaZOnYrc3NzKvk3Vi5Eao7i4mGlpabEOHTrIvc348eMZABYYGMhOnDjBNmzYwIyNjZmlpSXLyMjg6rm7uzNDQ0P22WefsQ0bNrCYmBgWEBDAALBt27Zx9XJycljLli2ZRCJhCxcuZCdPnmT79+9nU6ZMYWfPnmWMMVZSUsJ69+7NJBIJCw0NZTExMWzTpk2sUaNGzMHBgeXl5XHtWVtbs0aNGjErKyu2ZcsW9ttvv7ERI0YwAGzp0qWMMcays7NZVFQUA8Dmz5/PYmNjWWxsLEtOTuZid3d3l9l3X19fZm1tLVXm5+fHNm/ezGJiYlhMTAz7/vvvmVgsZqGhoVL1KmrzQ+7u7qxly5YsJyeHSSQStm7dOm5dhw4d2OjRo1lcXBwDwKKiorh158+fZxoaGszFxYXt2bOHHTp0iPXq1YsJBAK2e/durt6dO3eYlpYWc3BwYLt27WKHDx9mXl5ezMrKigFgiYmJXN2wsDAmEAjY2LFjWXR0NDtw4ABzc3NjEomE3blzh6u3YMECpuyPdmJiotTvrDwdO3ZkQqGQFRUVMcYYu3//Pvvmm2/Y7t272fnz51l0dDQbN24cU1NTY+fOneO2mz17NtPW1maFhYWMMcbS0tIYACYWi1lYWBhX75tvvmGmpqbcspeXFzM2NmaRkZHs/Pnz7NChQyw4OFjq/V28eDEDwIYPH86OHTvGfv75Z2Zra8v09PTYw4cPuXq+vr5MU1OT2dvbs4iICHb69GkWHBzMBAKB1N/OtGnTmFgsZtnZ2VL7vm7dOgaA3b59W+r9ev9vouxv/P3faXl/w/7+/kxDQ4NNnz6dnThxgu3cuZPZ2dkxU1NTlpaWVuH7/34fly9fZkVFRVKv4uJiqboTJkxgmpqaLC4ujjHG2JkzZ5iamhqbP3++VD15P1PW1tascePGzNHRke3atYv99ttvrEOHDkxDQ4MFBwezzp07swMHDrCDBw+y5s2bM1NTU6nvirLfgZWVFQsLC2OnTp1iISEhrEGDBqx///4yffn6+nLLubm5zMnJiRkZGbHly5ez06dPs1WrVjE9PT3m6enJSktLP/q+1QSUjGuQsi+hYcOGyVX/3r17DAALCAiQKr9y5QoDwObOncuVubu7MwDsypUrUnUdHByYl5cXt7xw4UIGgMXExFTY765duxgAtn//fqnysqT0fsKytrZmAoGA3bx5U6puz549ma6uLsvNzZXa9v0vr/djlzcZv6+kpIQVFRWxhQsXMkNDQ6kPpKLJuKw/V1dXxti7JAqAnT9/vtzYO3bsyExMTNjr16+5suLiYubo6MgaN27MxTJ06FAmFoulvmSLi4uZnZ2d1Bd3UlISa9CgAZs0aZJUfK9fv2ZmZmZsyJAhXBlfyXjo0KEMAHv+/Hm564uLi1lRURHr3r07GzRoEFd++vRpBoBduHCBMcbY9u3bmY6ODgsICGAeHh5cvc8++4z5+Phwy9ra2mzq1KkVxpOVlcXEYjHr27evVHlSUhITCoVSbfn6+jIA7Ndff5Wq27dvX9aiRQtu+fbt2wwAi4yMlKrXvn175uLiwi1XNhnHxsYyAGzZsmVS7ScnJzOxWMxmzZpV4f6+30d5L3V1dam6+fn5rG3btqxJkybs7t27zNTUlLm7u8sk7fd97DNlbW3NxGIxe/r0KVd28+ZNBoCZm5tzn3XGGDt06BADwI4cOSL1XgBgq1atkuozLCyMAWCXLl2S6uv9ZBweHs7U1NS4fyzK7Nu3jwFgv/3220fft5qAhqlrsXPnzgGAzKzC9u3bw97eHmfOnJEqNzMzQ/v27aXKWrdujSdPnnDLx48fR/PmzdGjR48K+42OjkbDhg3h7e2N4uJi7uXk5AQzMzOZGZstW7ZEmzZtpMp8fHyQk5ODGzduyLu7cjl79ix69OgBPT09qKurQ0NDA8HBwcjMzER6enqV2h47diyuXbuGhIQEbN68GU2bNsXnn38uUy83NxdXrlzBl19+CW1tba5cXV0do0aNwtOnT/HgwQMA736H3bt3h6mpqVS9oUOHSrV58uRJFBcXY/To0VLvuUgkgru7u8KzZBljUu0UFxcrtH1FbX5ow4YNcHZ2hkgkQoMGDaChoYEzZ85IDROXDdGXDS/HxMSgW7du6N27N/7880/k5eUhOTkZf//9t9TfZfv27bF161YsWrQIly9flhnaj42Nxdu3b2U+H5aWlvD09JT5fAgEAnh7e0uVffj5aNWqFVxcXBAVFcWV3bt3D1evXsXYsWPlfKcqFh0dDYFAgJEjR0r9bszMzNCmTRu5f88///wz4uLipF5XrlyRqiMUCvHrr78iMzMTzs7OYIxh165dUFdXl6qnyGfKyckJjRo14pbt7e0BvDvVpKWlJVP+/ntbZsSIEVLLPj4+AP73fVee6OhoODo6wsnJSep98/LyqjWzyCkZ1yBGRkbQ0tJCYmKiXPUzMzMBvJu5+SELCwtufRlDQ0OZekKhEG/fvuWWMzIy0Lhx44/2+/z5c7x69QqamprQ0NCQeqWlpeHFixdS9c3MzGTaKCv7MMaquHr1Knr16gXg3Yz0P/74A3FxcZg3bx4ASO1nZXz++ef47LPPsHHjRvzyyy8YO3ZsuZcQZWVlgTFW4e8F+N9+Z2ZmfvT9KfP8+XMA786vfvie79mzR+Y9/5Tff/9dpp3KXEb1vidPnkAoFMLAwAAAsHz5cnzzzTfo0KED9u/fj8uXLyMuLg69e/eW+l2IRCKpc71nzpxBz5490a1bN5SUlODixYvcOcz3k/GePXvg6+uLTZs2wc3NDQYGBhg9ejTS0tIAKP750NLSgkgkkioTCoXIz8+XKhs7dixiY2O5ORRRUVEQCoUYPny44m/aB54/fw7GGExNTWV+P5cvX5b792xvbw9XV1epl4uLi0y9Zs2aoWvXrsjPz8eIESNk3itFP1Nlv/sympqaHy3/8L1t0KCBzPeUPN8Vz58/x+3bt2XeMx0dHTDGFP588IFmU9cg6urq6N69O44fP46nT59+MimW/dGmpqbK1E1JSYGRkZHCMRgbG8tM/vqQkZERDA0NK5zVraOjI7Vc9uVYXll5/yB8SCQSITs7W6b8ww/Y7t27oaGhgejoaKkv1UOHDn2yD3mNGTMG8+fPh0AggK+vb7l19PX1oaamhtTUVJl1KSkpAMD9bgwNDT/6/pQpq79v3z5YW1tXaR8AwMXFRWpyFPC/fxQq49mzZ7h+/Trc3d3RoMG7r5Xt27ejW7duWL9+vVTd169fy2zfvXt3BAcH4+rVq3j69Cl69uwJHR0dtGvXDjExMUhJSUHz5s1haWnJbWNkZISVK1di5cqVSEpKwpEjRzB79mykp6fjxIkTUp+PD1X28wEAw4cPR1BQELZu3YqwsDD88ssvGDhwIPT19SvV3vuMjIwgEAhw8eLFcidwlldWFZs2bcKxY8fQvn17rF27FkOHDkWHDh249dXxmXpfcXExMjMzpb4X5PmuMDIyglgsxpYtWypcX9PRkXENM2fOHDDG4O/vj8LCQpn1RUVFOHr0KADA09MTwLsvvffFxcXh3r17MtfEyqNPnz54+PAhzp49W2Gd/v37IzMzEyUlJTL/fbu6uqJFixZS9e/cuYNbt25Jle3cuRM6OjpwdnYG8L8vmfKOXm1sbPDw4UNuVivw7r/kP//8U6qeQCBAgwYNpIbZ3r59i19++UXOvf80X19feHt7Y+bMmVLDce+TSCTo0KEDDhw4ILU/paWl2L59Oxo3bozmzZsDADw8PHDmzBnuyBcASkpKsGfPHqk2vby80KBBAzx69Kjc91zRy1h0dHRkti87WlHU27dv8dVXX6G4uBizZs3iygUCgUzyuH37NmJjY2Xa6NGjB4qLi/Hdd9+hcePG3DWyPXr0wOnTp7mh0opYWVkhMDAQPXv25E59uLm5QSwWy3w+nj59Wu414/LS19fHwIED8fPPPyM6OhppaWlKGaIG3n22GGN49uxZub/jVq1aKaUfAEhISMDkyZMxevRoXLx4Ea1bt8bQoUORlZXF1amOz9SHduzYIbW8c+dOAPjoTXr69++PR48ewdDQsNz3TZk3VlEVOjKuYdzc3LB+/XoEBATAxcUF33zzDVq2bImioiLEx8cjMjISjo6O8Pb2RosWLTB+/HisWbMGampq6NOnDx4/fozvvvsOlpaWmDZtmsL9T506FXv27MGAAQMwe/ZstG/fHm/fvsXvv/+O/v37w8PDA8OGDcOOHTvQt29fTJkyBe3bt4eGhgaePn2Kc+fOYcCAARg0aBDXpoWFBf7zn/8gJCQE5ubm2L59O2JiYvDjjz9y55GaNm0KsViMHTt2wN7eHtra2rCwsICFhQVGjRqFjRs3YuTIkfD390dmZiaWLFkCXV1dqdj79euH5cuXw8fHB+PHj0dmZiYiIiKUejRhYWEh11FBeHg4evbsCQ8PD8yYMQOamppYt24d/vrrL+zatYsb3p4/fz6OHDkCT09PBAcHQ0tLC//9739lLsewsbHBwoULMW/ePPz777/o3bs39PX18fz5c1y9ehUSiQShoaFK28+KJCUl4fLlyygtLUV2djZ3048nT55g2bJl3JAm8O4L8vvvv8eCBQvg7u6OBw8eYOHChWjSpInMOWoXFxfo6+vj1KlTGDNmDFfeo0cPfP/999zPZbKzs+Hh4QEfHx/Y2dlBR0cHcXFxOHHiBL744gsAQMOGDfHdd99h7ty5GD16NIYPH47MzEyEhoZCJBJhwYIFlX4fxo4diz179iAwMBCNGzf+6D8KiujcuTPGjx+PMWPG4Nq1a/j8888hkUiQmpqKS5cuoVWrVvjmm28+2c5ff/1V7jyApk2bwtjYGLm5uRgyZAiaNGmCdevWQVNTE7/++iucnZ0xZswY7m+8Oj5T79PU1MSyZcvw5s0btGvXDn/++ScWLVqEPn36oEuXLhVuN3XqVOzfvx+ff/45pk2bhtatW6O0tBRJSUk4deoUpk+fLnXEXyPxN3eMfMzNmzeZr68vs7KyYpqamkwikbC2bduy4OBglp6eztUrKSlhP/74I2vevDnT0NBgRkZGbOTIkdxlQWXenxX8vvJmJGdlZbEpU6YwKysrpqGhwUxMTFi/fv3Y/fv3uTpFRUUsIiKCtWnTholEIqatrc3s7OzY119/zf7++2+unrW1NevXrx/bt28fa9myJdPU1GQ2NjZs+fLlMrHs2rWL2dnZMQ0NDQaALViwgFu3bds2Zm9vz0QiEXNwcGB79uwpN/YtW7awFi1aMKFQyGxtbVl4eDjbvHmzzCzWysymrkhFM8EvXrzIPD09mUQiYWKxmHXs2JEdPXpUZvs//viDuyzIzMyMzZw5k0VGRsrEzNi7WageHh5MV1eXCYVCZm1tzb788kt2+vRpro4qZ1PjvZm5+vr6zMXFhU2dOlXq0qoyBQUFbMaMGaxRo0ZMJBIxZ2dndujQoQpnwQ8aNIgBYDt27ODKCgsLmUQiYWpqaiwrK4srz8/PZxMmTGCtW7dmurq6TCwWsxYtWrAFCxZIzdpljLFNmzax1q1bM01NTaanp8cGDBggE6+vry+TSCQyMVX0XpaUlDBLS0sGgM2bN6/C96sylzYx9u7vuEOHDtzfTtOmTdno0aPZtWvXZOq+72OzqQGwn376iTHG2MiRI5mWlpbM+7B3714GgK1YsUIqFnk+U2Wf9Q8BYBMnTiz3/Xl/dn7Z7+D27dusW7duTCwWMwMDA/bNN9+wN2/eSG3/4Wxqxhh78+YNmz9/PmvRogX3u27VqhWbNm3aJy8JqwkEjFXxbgiEfISNjQ0cHR25C/8JIaQ8fn5+2LdvH968ecN3KLygc8aEEEIIzygZE0IIITyjYWpCCCGEZ3RkTAghhPCMkjEhhBDCM0rGhBBCCM/oph+VUFpaipSUFOjo6JR7b2JCCCF1H2MMr1+/hoWFBdTUqnZsS8m4ElJSUqTukUsIIaT+Sk5O/uSzBD6FknEllD0I4dq1ZGhr636iNiHKl5BwHTt2bMLVqxeRkfEcWloSODm1g7//NLi6dlJZv2/evMb69Utw714C7t27jaysTEycOBuTJs0pt35u7husWrUIx48fRHZ2Fmxtm8Pffyr69ftS6X1dvx6LjRuX4ebNqygoKICZmQUGDBiOgIBZ5dYvc+DADsydG4C9e8+hVSvnT78JhPx/b97kwNXVUubhOJVBybgSyoamtbV1oaNDyZhUryVLvsOaNYvRqZMHZs78Ho0aWeHZsyRERi7H6NH9sGrVz/jiixGfbqgSXr16ib17t8HBoQ369BmEnTs3QSgUVvg5GD/+S9y6FYc5c36ArW1zHDq0E9Onj4NQKMKgQT5K6+vgwZ2YPHkUvL2HYNWqXyCRaOPJk0d4/jzlk59RkUgMAJBItOnzTCpFGacrKRkTUotERCzAqlWL8N13SzFhwgypdQMGDIOnZ0vMmzcRPXr0h66untL7b9zYGnfvZkEgEODlyxfYuXNThXXPnPkNFy7E4L//3YmBA98967dzZw88ffoEixbNxH/+M1TmQfaV6Ss19RlmzRqPkSO/Rnj4Oq68c2ePSu4lIdWPZlMTUktcv34Zq1YtwuDBvjKJGHj3GMoRI8YjJycbly6dUUkMAoFA7qOAEycOQiLRRv/+g6XKhw4dg7S0FNy4cUUpfe3atQl5ebmYOPFbueIipCaiZExILbFy5fcQCASYOXNhhXWsrGwBAKmpT2XWMcZQXFws10sZ7t//C599Zo8GDaQH4OztWwMAHjz4Syn9XL58AQ0bGuCff+6jZ08nWFk1QOvWJvj22wl4/TpHKX0Qomo0TE1ILZCTk40LF07B3b0XGjWyqrBebu67J95oaUlk1sXG/o7Bg+Ubur18ORGWljaVirVMVlYmrK1tZcobNjTg1itDWtoz5Ofn4euvByMwcA5cXFbi1q04REQswIMHf+HgwYt0CSKp8SgZE1IL3L+fgOLiYtjZtfpovevXYwH87+jzfa1bu+C33+Lk6s/U1ELxIMvxsSSorARZWlqK/Px8zJmzAIGBswEAnTp1g4aGJhYsmIqLF8/g8897KKUvQlSFkjEhtUDZcKuhoXGFdd68eY3Dh3fByqoJ2rRxlVkvkWijZUsnufr7cGi5MvT1Dcs9+n316iWA/x0hK6OfxMS/0a2bl1S5h0cfLFgwFX/9dYOSManxKBmTeuvZsyS8fPmCl74NDIw+Otz8IXPzdzcUSE5+XGGddeuW4M2b1/j++zXlHnVW9zC1vX0rHDq0C8XFxVLJ/f79BABAixaOVWr/f/20xo0bl2XKyx5IV9U7IxFSHSgZk3rp2bMkuLvb4+3bPF76F4u18Pvv9+ROyPb2rWBj0xSHDu3CzJnfo2FDfan1+/dvx5o1i+HtPQRDhviW20Z1D1P37j0IO3b8hGPH9mPAgKFc+d6922BmZgFn5w5V7gMA+vX7P+zYEYmzZ4/D0bEtV3727G8AAGfnjkrphxBVomRM6qWXL1/g7ds8hISEwMbGplr7fvz4MUJCQvDy5Qu5k7FAIMCSJT9h9Oi+6NevHQICvoWNTTNkZKTh8OHdOHXqCIYM8cOSJZEVtqGtrVPu8LWizp49jry8XOTmvgYAPHx4F9HR+wAA3bv3hVisBQDw9OyDzz/viblzv8GbNzmwsWmGw4d34dy5E1izZjt3jXFs7O8YOrQ7pk0LxrRpwQr35e7eCz17emPlyoUoLS2Fs3NH3L59DStWhKJHj/5o376LXPv1xx9nyx15eH+fCFEVSsakXrOxsYGdnR3fYcilc2cPHD16BatXh2Hp0u+QmZmB0tJSNG5sjV27YqrtvOicOd/g6dMn3HJ09F5ER+8FIDu8vWnTAfz44zxERATj1auXaNrUDuvW7cKAAcO4OowxlJSUoLS0tNJ9rV+/BytWhGLHjkisWBEKU1MLfPXVNAQFLZB7v8LCyr9OWRlD9oR8ioCVnVipZ9atW4elS5ciNTUVLVu2xMqVK9G1a1e5ts3JyYGenh7u38+m2+fVUgkJN9C7twu2bt1a7cn4/v378PPzw4kT16t8L+RJk0bh6NE9OHz4T6Uc9RJC5Pf6dQ7s7PSQnZ0NXd2q5YJ6ObNhz549mDp1KubNm4f4+Hh07doVffr0QVJSEt+hEaKQxYv/CzOzRggMHIG8vFy+wyGEVFK9HKZevnw5xo0bh6+++goAsHLlSpw8eRLr169HeHi43O3k5eV+9N66pObKz38LACgoKMDbt2+rte+CggIuhqomUHV1dZw9+787WVFCJqT6KPPzVu+GqQsLC6GlpYW9e/di0KBBXPmUKVNw8+ZN/P777zLbFBQUcF+gwLthanqeMSGEEAA0TF0ZL168QElJCUxNTaXKTU1NkZaWVu424eHh0NPT416UiAkhhChTvRymBmRvxccYq/D2fHPmzEFQUBC3XHZk/PXXX0MoFKo0TqI6OTk5yM/PBwB4eFTP4/bOnTsHABCJRFX+T5oQwq+CggJs3LhRKW3Vu2RsZGQEdXV1maPg9PR0maPlMkKhsNykq6GhAQ0NDZXESVTP0NCQ+7lNmzbV0ufDhw+rpR9CiOqVdzleZdW7YWpNTU24uLggJiZGqjwmJgadOnXiKSpCCCH1Wb07MgaAoKAgjBo1Cq6urnBzc0NkZCSSkpIwYcIEvkMjhBBSD9XLZDx06FBkZmZi4cKFSE1NhaOjI3777TdYW1vzHRohhJB6qF4mYwAICAhAQEAA32EQQggh9e+cMSGEEFLTUDImhBBCeEbJmBBCCOEZJWNCCCGEZ5SMCSGEEJ5RMiaEEEJ4RsmYEEII4RklY0IIIYRnlIwJIYQQnlEyJoQQQnhGyZgQQgjhGSVjQgghhGeUjAkhhBCeUTImhBBCeEbJmBBCCOEZJWNCCCGEZ5SMCSGEEJ5RMiaEEEJ4RsmYEEII4RklY0IIIYRnlIwJIYQQnlEyJoQQQnhGyZgQQgjhGSVjQgghhGd1Jhk/fvwY48aNQ5MmTSAWi9G0aVMsWLAAhYWFUvWSkpLg7e0NiUQCIyMjTJ48WaYOIYQQUp0a8B2Asty/fx+lpaXYuHEjmjVrhr/++gv+/v7Izc1FREQEAKCkpAT9+vWDsbExLl26hMzMTPj6+oIxhjVr1vC8B4QQQuqrOpOMe/fujd69e3PLtra2ePDgAdavX88l41OnTuHu3btITk6GhYUFAGDZsmXw8/NDWFgYdHV1eYmdEEJI/SZXMg4KClK44fnz58PAwEDh7ZQpOztbKobY2Fg4OjpyiRgAvLy8UFBQgOvXr8PDw6PcdgoKClBQUMAt5+TkqC5oQggh9Y5cyXjlypVwc3ODpqamXI1eunQJgYGBvCbjR48eYc2aNVi2bBlXlpaWBlNTU6l6+vr60NTURFpaWoVthYeHIzQ0VGWxEkIIqd/kHqY+ePAgTExM5Kqro6NT6YA+FBIS8slEGBcXB1dXV245JSUFvXv3xuDBg/HVV19J1RUIBDLbM8bKLS8zZ84cqdGBnJwcWFpayrsLhBBCyEfJlYyjoqKgp6cnd6MbN26UOQKtrMDAQAwbNuyjdWxsbLifU1JS4OHhATc3N0RGRkrVMzMzw5UrV6TKsrKyUFRU9NF4hUIhhEKh4sETQgghcpArGfv6+irUqI+PT6WCKY+RkRGMjIzkqvvs2TN4eHjAxcUFUVFRUFOTvnLLzc0NYWFhSE1Nhbm5OYB3k7qEQiFcXFyUFjMhhBCiiCrNpn7z5g1KS0ulyviakZySkoJu3brBysoKERERyMjI4NaZmZkBAHr16gUHBweMGjUKS5cuxcuXLzFjxgz4+/vTTGpCCCG8UTgZJyYmIjAwEOfPn0d+fj5XXnbetaSkRKkByuvUqVP4559/8M8//6Bx48ZS6xhjAAB1dXUcO3YMAQEB6Ny5M8RiMXx8fLhLnwghhBA+CFhZppJTp06dAABTpkyBqampzMQnd3d35UVXQ+Xk5EBPTw+BgYF0LrkO2L59u9QIT+fOnTF79mypOn5+fsjMzPxkW4GBgfDy8uKWnzx5gsDAQG757du3kEgkGDlypBIiJ4TwqaCgAGvXrkV2dnaVR1cVPjK+ffs2rl+/jhYtWlSp47pgVHdTaEtEfIdBqmj7tnxkvMjmlsu7jjwzM1Pq1EdF3h8tAt7d9e3D7SQidYz1MqtktISQmuJNbj7WrlVOWwon43bt2iE5OZmSMalz1NTUYGhoWO5/uIaGhnK1IRJJ/3Omrq4OY2NjAO8S+odzLAghBKhEMt60aRMmTJiAZ8+ewdHRERoaGlLrW7durbTgCKlOhoaGOHr0aLnrtm7dWqk2ra2tuTa9vb3lOromhNQ/CifjjIwMPHr0CGPGjOHKBAIB7xO4CCGEkNpK4WQ8duxYtG3bFrt27Sp3AhchhBBCFKNwMn7y5AmOHDmCZs2aqSIeQgghpN5R+3QVaZ6enrh165YqYiGEEELqJYWPjL29vTFt2jQkJCSgVatWMhO4/vOf/ygtOEKqQ99uTnhRoKvyu7BFRUVB/W0SxMWpKu2HEFL7KHzTjw/v9yzVWD2ZwFV2048rB7+n64zriNe6XaulH82CJxAWJHHLB0/FYX7Er9izdjIcm6vmSWAxlxJw8sJt/PUwGekvsmGor4O2DjaYOLonrBsZq6RPQuqDN7n56DDoO35u+kHXSRJSu2zecw5GBjr4enh3NDY3QFr6K0TuPosvA1Zh16pANLOhG5AQwrcqPSiCEFLz/XfhWBjqa0uVdWjbDD1HhePnAxexMGgwT5ERQsrINYFr9erVMrf5+5gNGzbg9evXlQ6KkLro0KFD2L4nGr8euyz3NldvPULLXjNx7Fw8VkUdR7dh36P9wPkY9+1GJCany9XGh4kYAEwM9WBmpIfUjFdyx0IIUR25kvG0adMUSq6zZs2iOw2RWqP/2CXw9PTE0KFDVdrP5s2bsXzdL9iw47TC267achwpz7OwMOhLhEz5Ek+evcDE4CiUlFTutFFyaiZS0rPQzNq0UtsTQpRLrmFqxhi6d++OBg3kG9V++/ZtlYIipDrl5RciLy8PeXl5fIdSoabWpvhxtg+3rK4uQNCi7fjrYTLa2Fsr1FZxSQm+W74XWiIhRn/xubJDJYRUglzZdcGCBQo1OmDAABgYGFQqIEKILI+ODlLLzZuYAwBSnmehjb01GGMo+WByZQN1dZl2GGP4btle3EhIxMrgUTA3aaiymAkh8lNJMiaEKJeerkRqWVPj3Uc3v7AIAHAo5hrmR/wqVefOqaVSy4wxBC/fi+izN7B4xlB4dnJUYcSEEEXQbGpC6gCPjg7Ys3ZyhevLEvHBU9fwfdBgePdwqcboCCGfQsmYkDqgoa4EDT84ei7DGMOCFftw8NQ1LJjyfxjk1a6aoyOEfAolY0JqgCs3/8GztCyZcrFIo5zailm87jD2n7iKL7zaoXkTM9y694Rbp6nRAPbNGlW5D0JI1VAyJqQGWL7pt3LLF80YUuW2z1++CwA4cDIOB07GSa2zMNVHzC9zq9wHIaRqFE7GCxcuxIwZM6ClpSVV/vbtWyxduhTBwcFKC46Qum5Qr3YY1Ovjw8blrW9kZiAzQasilGwJqfkUfoRiaGgo3rx5I1Oel5eH0NBQpQRFSF1kZWUFW5vGsG5kxHcohJAaRuEjY8YYBAKBTPmtW7fo2mJSKwVP/gKv1JtBKBSqtJ///ve/Mk9tIoQQQIFkrK+vD4FAAIFAgObNm0sl5JKSErx58wYTJkxQSZCEqFK3jg7V9ghFQggpj9zJeOXKlWCMYezYsQgNDYWenh63TlNTEzY2NnBzc1NJkIoqKChAhw4dcOvWLcTHx8PJyYlbl5SUhIkTJ+Ls2bMQi8Xw8fFBREQENDU1+QuYEEJIvSZ3Mvb19QUANGnSBJ06dYKGRtUvuVCVWbNmwcLCArdu3ZIqLykpQb9+/WBsbIxLly4hMzMTvr6+YIxhzZo1PEVLCCGkvlP4nLG7uztKS0vx8OFDpKeno/SD++F+/jm/N54/fvw4Tp06hf379+P48eNS606dOoW7d+8iOTkZFhYWAIBly5bBz88PYWFh0NXV5SNkwrM7D58iW5gADQ0N2NnZqayf4OBg5GSlwVBbHUvm+Hx6A0JIvaFwMr58+TJ8fHzw5MkTMMak1gkEApSUlCgtOEU9f/4c/v7+OHTokMylVwAQGxsLR0dHLhEDgJeXFwoKCnD9+nV4eHiU225BQQEKCgq45ZycHOUHT3gzKWQrnr/IhrGxMY4ePaqyfuLj45GRkQFTI71PVyaE1CsKX9o0YcIEuLq64q+//sLLly+RlZXFvV6+fKmKGOXCGIOfnx8XX3nS0tJgair9/FZ9fX1oamoiLS2twrbDw8Ohp6fHvSwtLZUaOyGEkPpN4WT8999/Y/HixbC3t0fDhg2lktT7k7qUJSQkhJvFXdHr2rVrWLNmDXJycjBnzpyPtlfeZVkVXa5VZs6cOcjOzuZeycnJVd4vQgghpIzCw9QdOnTAP//8g2bNmqkiHhmBgYEYNmzYR+vY2Nhg0aJFuHz5ssy1oq6urhgxYgS2bdsGMzMzXLlyRWp9VlYWioqKZI6Y3ycUClV+DSohhJD6S65kfPv2be7nSZMmYfr06UhLS0OrVq1kZlW3bt1aqQEaGRnByOjTdyxavXo1Fi1axC2npKTAy8sLe/bsQYcOHQAAbm5uCAsLQ2pqKszN3z2c/dSpUxAKhXBxoUfKEUII4YdcydjJyQkCgUBqwtbYsWO5n8vW8TmBy8rKSmpZW1sbANC0aVM0btwYANCrVy84ODhg1KhRWLp0KV6+fIkZM2bA39+fZlITQgjhjVzJODExUdVxVAt1dXUcO3YMAQEB6Ny5s9RNPwghhBC+yJWMra2tVR2H0tnY2MhcegW8O4KOjo7mISJCCCGkfApP4Dpy5Ei55QKBACKRCM2aNUOTJk2qHBghhBBSXyicjAcOHChz/hiQPm/cpUsXHDp0CPr6+koLlJDabsCAAXibnYqGokK+QyGE1DAKX2ccExODdu3aISYmhrvuNiYmBu3bt0d0dDQuXLiAzMxMzJgxQxXxEqJ0RzfNwJkzZ7B7926V9vPVV19heuBoBIzqpdJ+CCG1j8JHxlOmTEFkZCQ6derElXXv3h0ikQjjx4/HnTt3sHLlSqnZ1oTUZBItEUolEr7DIITUYwofGT969Kjcy4B0dXXx77//AgA+++wzvHjxourREUIIIfWAwsnYxcUFM2fOREZGBleWkZGBWbNmoV27dgDe3TKz7NpeQgghhHycwsPUmzdvxoABA9C4cWNYWlpCIBAgKSkJtra2OHz4MADgzZs3+O6775QeLCGqsHXf73hZch8SiQQ+Pqp7tKG3tzf31KazO+errB9CSO2jcDJu0aIF7t27h5MnT+Lhw4dgjMHOzg49e/aEmtq7A+2BAwcqO05CVObnAxe5RyiqMhkTQkhFFE7GwLvLmHr37o3evXsrOx5CCCGk3pErGa9evRrjx4+HSCTC6tWrP1p38uTJSgmMEEIIqS/kSsYrVqzAiBEjIBKJsGLFigrrCQQCSsaEEEKIghR+UERdeWgEIYQQUlMofGlTmcLCQjx48ADFxcXKjIcQQgipdxROxnl5eRg3bhy0tLTQsmVLJCUlAXh3rviHH35QeoCEEEJIXadwMp4zZw5u3bqF8+fPQyQSceU9evTAnj17lBocIYQQUh8ofGnToUOHsGfPHnTs2BECgYArd3BwwKNHj5QaHCGEEFIfKJyMMzIyYGJiIlOem5srlZwJqS3smzWCsZklGjZsqNJ+QkJCwPKeQoIslfZDCKl9FE7G7dq1w7FjxzBp0iQA4BLwTz/9BDc3N+VGR0g1+O/CMXit21Xl/bi4uECzwAjCgiSV90UIqV0UTsbh4eHo3bs37t69i+LiYqxatQp37txBbGwsfv/9d1XESAghhNRpCk/g6tSpE/744w/k5eWhadOmOHXqFExNTREbGwsXFxdVxEgIIYTUaZW6N3WrVq2wbds2ZcdCSJ12/fp17pxx+zZN+Q6HEFKDVCoZl5aW4p9//kF6ejpKS0ul1n3++edKCYyQ6jIxOAov3mxDw4YNERERobJ+QkJC6BGKhJByKZyML1++DB8fHzx58gSMMal1AoEAJSUlSguOkOpw759n3CMUCSGEDwon4wkTJsDV1RXHjh2Dubk5Xc5ECCGEVJHCE7j+/vtvLF68GPb29mjYsCH09PSkXnw7duwYOnToALFYDCMjI3zxxRdS65OSkuDt7Q2JRAIjIyNMnjwZhYWFPEVLCCGEVOLIuEOHDvjnn3/QrFkzVcRTJfv374e/vz8WL14MT09PMMaQkJDArS8pKUG/fv1gbGyMS5cuITMzE76+vmCMYc2aNTxGTgghpD6TKxnfvn2b+3nSpEmYPn060tLS0KpVK2hoaEjVbd26tXIjlFNxcTGmTJmCpUuXYty4cVx5ixYtuJ9PnTqFu3fvIjk5GRYWFgCAZcuWwc/PD2FhYdDV1a32uAkhhBC5krGTkxMEAoHUhK2xY8dyP5et43MC140bN/Ds2TOoqamhbdu2SEtLg5OTEyIiItCyZUsAQGxsLBwdHblEDABeXl4oKCjA9evX4eHhUW7bBQUFKCgo4JZzcnJUuzOEEELqFbmScWJioqrjqLJ///0XwLvLR5YvXw4bGxssW7YM7u7uePjwIQwMDJCWlgZTU1Op7fT19aGpqYm0tLQK2w4PD0doaKhK4yeEEFJ/yZWMra2tVR1HhUJCQj6ZCOPi4rjrnefNm4f/+7//AwBERUWhcePG2Lt3L77++msAKHf2d9lRfUXmzJmDoKAgbjknJweWlpYK7wshhBBSnkrd9KM6BQYGYtiwYR+tY2Njg9evXwN49yjHMkKhELa2tkhKendjfjMzM1y5ckVq26ysLBQVFckcMb9PKBRCKBRWdhcIIYSQj6rxydjIyAhGRkafrOfi4gKhUIgHDx6gS5cuAICioiI8fvyYO7J3c3NDWFgYUlNTYW5uDuDdpC6hUEj31a7HRn/RFS9LTCGRSFTaz9GjR6FZ8ISe2kQIkVHjk7G8dHV1MWHCBCxYsACWlpawtrbG0qVLAQCDBw8GAPTq1QsODg4YNWoUli5dipcvX2LGjBnw9/enmdT1mN+X7tXyCEVCCKlInUnGALB06VI0aNAAo0aNwtu3b9GhQwecPXsW+vr6AAB1dXUcO3YMAQEB6Ny5M8RiMXx8fFR6P2JCCCHkUwTswxtMy+HVq1fYt28fHj16hJkzZ8LAwAA3btyAqakpGjVqpIo4a5ScnBzo6enhysHvoS0R8R0OUYLqOjKmYWpC6o43ufnoMOg7ZGdnV3l0VeEj49u3b6NHjx7Q09PD48eP4e/vDwMDAxw8eBBPnjzBzz//XKWACKluuXn5yFXPBQCVnjfetGkT3manoqGoEAGjeqmsH0JI7aPwvamDgoLg5+eHv//+GyLR/44K+/TpgwsXLig1OEKqg/dXEejevfsnZ+1X1eHDh7Fj7zHsO35Vpf0QQmofhZNxXFwcd83u+xo1avTRG2cQQgghpHwKJ2ORSFTu7SAfPHhAz4MlhBBCKkHhZDxgwAAsXLgQRUVFAN7d0SopKQmzZ8/m7nxFCCGEEPkpnIwjIiKQkZEBExMTvH37Fu7u7mjWrBl0dHQQFhamihgJIYSQOk3h2dS6urq4dOkSzp49ixs3bqC0tBTOzs7o0aOHKuIjhBBC6jyFk/Hjx49hY2MDT09PeHp6qiImQgghpF5ReJja1tYWXbp0wcaNG/Hy5UtVxEQIIYTUKwon42vXrsHNzQ2LFi2ChYUFBgwYgL1796KgoEAV8RFCCCF1nsLJ2NnZGUuXLkVSUhKOHz8OExMTfP311zAxMcHYsWNVESMhdULbtm3RsV1ruLay5TsUQkgNU6l7U3/oxo0bGDduHG7fvo2SkhJlxFWj0b2p65Y7D58iW+gIDQ0N2NnZqbQvujc1IXUHr/emLpOcnIxdu3Zh586dSEhIgJubG9auXVulYAjhQ8vmjfFatxXfYRBC6jGFk3FkZCR27NiBP/74Ay1atMCIESNw6NAh2NjYqCA8QgghpO5TOBl///33GDZsGFatWgUnJycVhEQIIYTULwon46SkJAgEAlXEQggvzl++i1fqhRAKhejSpYvK+pk4cSKyMp/DSE+IqKUTVNYPIaT2kSsZ3759G46OjlBTU0NCQsJH67Zu3VopgRFSXRauPoDnL7JhbGys0mSclJSEjIwM5BrpqawPQkjtJFcydnJyQlpaGkxMTODk5ASBQID3J2GXLQsEgnoxm5oQQghRJrmScWJiIvd4xMTERJUGRAghhNQ3ciVja2tr7ucnT56gU6dOaNBAetPi4mL8+eefUnUJIYQQ8mkK34HLw8Oj3HtSZ2dnw8PDQylBEUIIIfWJwsm47NzwhzIzMyGRSJQSFCGEEFKfyH1p0xdffAHg3WQtPz8/CIVCbl1JSQlu376NTp06KT9CQgghpI6TOxnr6b27HIMxBh0dHYjFYm6dpqYmOnbsCH9/f+VHSAghhNRxcifjqKgoAICNjQ1mzJhRI4ekHz58iJkzZ+KPP/5AYWEhWrVqhUWLFkmdy05KSsLEiRNx9uxZiMVi+Pj4ICIiApqamjxGTgghpD5T+A5cCxYsUEUcStGvXz80b96cS7QrV65E//798ejRI5iZmaGkpAT9+vWDsbExLl26hMzMTPj6+oIxhjVr1vAdPuGJlkgTWlpa0NLSUmk/48aNQ2HOM+ho5Km0H0JI7VOpRyju27cPv/76K5KSklBYWCi17saNG0oLThEvXryAsbExLly4gK5duwIAXr9+DV1dXZw+fRrdu3fH8ePH0b9/fyQnJ8PCwgIAsHv3bvj5+SE9PV3uR2DRIxTrnte6XaulH3qEIiF1hzIfoajwbOrVq1djzJgxMDExQXx8PNq3bw9DQ0P8+++/6NOnT5WCqQpDQ0PY29vj559/Rm5uLoqLi7Fx40aYmprCxcUFABAbGwtHR0cuEQOAl5cXCgoKcP369QrbLigoQE5OjtSLEEIIURaFk/G6desQGRmJtWvXQlNTE7NmzUJMTAwmT56M7OxsVcQoF4FAgJiYGMTHx0NHRwcikQgrVqzAiRMn0LBhQwBAWloaTE1NpbbT19eHpqYm0tLSKmw7PDwcenp63MvS0lKVu0IIIaSeUTgZJyUlcZcwicVivH79GgAwatQo7Nq1S7nRAQgJCYFAIPjo69q1a2CMISAgACYmJrh48SKuXr2KAQMGoH///khNTeXaK+8a6YqunS4zZ84cZGdnc6/k5GSl7yep+168eIHn6ZnIyKSRFUKINIUncJmZmSEzMxPW1tawtrbG5cuX0aZNGyQmJqISp58/KTAwEMOGDftoHRsbG5w9exbR0dHIysrixu7XrVuHmJgYbNu2DbNnz4aZmRmuXLkitW1WVhaKiopkjpjfJxQKpa6rJnVLRGQ0XhRcgK6uLiZNmqSyfsaMGYOMjAyYGunh7M75KuuHEFL7KJyMPT09cfToUTg7O2PcuHGYNm0a9u3bh2vXrnE3BlEmIyMjGBkZfbJeXt67GapqatIH+2pqaigtLQUAuLm5ISwsDKmpqTA3NwcAnDp1CkKhkDuvTOqf387f5B6hqMpkTAghFVE4GUdGRnLJbcKECTAwMMClS5fg7e2NCRP4e2C6m5sb9PX14evri+DgYIjFYvz0009ITExEv379AAC9evWCg4MDRo0ahaVLl+Lly5eYMWMG/P39qzwTjhBCCKkshZOxmpqa1NHnkCFDMGTIEKUGVRlGRkY4ceIE5s2bB09PTxQVFaFly5Y4fPgw2rRpAwBQV1fHsWPHEBAQgM6dO0vd9IMQQgjhi1zJ+Pbt23I32Lp160oHU1Wurq44efLkR+tYWVkhOjq6miIihBBCPk2uZOzk5ASBQPDJCVoCgQAlJSVKCYwQQgipL+RKxomJiaqOgxBCCKm35ErG1tbWqo6DEEIIqbcUvukHAPzyyy/o3LkzLCws8OTJEwDAypUrcfjwYaUGRwghhNQHCifj9evXIygoCH379sWrV6+4c8QNGzbEypUrlR0fIYQQUucpfGnTmjVr8NNPP2HgwIH44YcfuHJXV1fMmDFDqcERUh0+b2+HzHyJyq81X7t2LdTeJkNc8lyl/RBCah+Fk3FiYiLatm0rUy4UCpGbm6uUoAipTiFTv6yWRyhaW1tDswAQFij/trGEkNpN4WHqJk2a4ObNmzLlx48fh4ODgzJiIoQQQuoVhY+MZ86ciYkTJyI/Px+MMVy9ehW7du1CeHg4Nm3apIoYCSGEkDpN4WQ8ZswYFBcXY9asWcjLy4OPjw8aNWqEVatWffLpSoTUZydPnkTxm2fQVn+D/p6yp3oIIfWXwskYAPz9/eHv748XL16gtLQUJiYmAIBnz56hUaNGSg2QEFUbMnEV0l8tgaGhIbZu3aqyftauXcs9QpGSMSHkfZW6zriMkZERTExMkJaWhkmTJqFZs2bKiouQavMi6zUyMjKQmZnJdyiEkHpK7mT86tUrjBgxAsbGxrCwsMDq1atRWlqK4OBg2Nra4vLly9iyZYsqYyWEEELqJLmHqefOnYsLFy7A19cXJ06cwLRp03DixAnk5+fj+PHjcHd3V2WchBBCSJ0ldzI+duwYoqKi0KNHDwQEBKBZs2Zo3rw53XWLEEIIqSK5h6lTUlK464htbW0hEonw1VdfqSwwQgghpL6QOxmXlpZCQ0ODW1ZXV4dEIlFJUIQQQkh9IvcwNWMMfn5+EAqFAID8/HxMmDBBJiEfOHBAuRESQgghdZzcydjX11dqeeTIkUoPhhBCCKmP5E7GUVFRqoyDkDrP0NAQApTAqKEW36EQQmqYSt2Bi5C6JOirfnglsIZIJFJpP1u3boVmwRMIC5JU2g8hpPahZEzqvf6ebavlEYqEEFKRKt0OkxBCCCFVR8mYEEII4VmtScZhYWHo1KkTtLS00LBhw3LrJCUlwdvbGxKJBEZGRpg8eTIKCwul6iQkJMDd3R1isRiNGjXCwoULwRirhj0gNVVicjr+/fdfPHnyRKX9/PDDD5gVvBwhK/eptB9CSO1Ta84ZFxYWYvDgwXBzc8PmzZtl1peUlKBfv34wNjbGpUuXkJmZCV9fXzDGsGbNGgBATk4OevbsCQ8PD8TFxeHhw4fw8/ODRCLB9OnTq3uXSA0x7ttIPH+RDWNjYxw9elRl/fzxxx/cIxQJIeR9tSYZh4aGAkCFz5s9deoU7t69i+TkZFhYWAAAli1bBj8/P4SFhUFXVxc7duxAfn4+tm7dCqFQCEdHRzx8+BDLly9HUFAQBAJBde0OIYQQwqk1w9SfEhsbC0dHRy4RA4CXlxcKCgpw/fp1ro67uzt3F7GyOikpKXj8+HGFbRcUFCAnJ0fqRQghhChLnUnGaWlpMDU1lSrT19eHpqYm0tLSKqxTtlxWpzzh4eHQ09PjXpaWlkqOnhBCSH3GazIOCQmBQCD46OvatWtyt1feMDNjTKr8wzplk7c+NkQ9Z84cZGdnc6/k5GS5YyKEEEI+hddzxoGBgRg2bNhH69jY2MjVlpmZGa5cuSJVlpWVhaKiIu7o18zMTOYIOD09HQBkjpjfJxQKpYa2CSGEEGXiNRkbGRnByMhIKW25ubkhLCwMqampMDc3B/BuUpdQKISLiwtXZ+7cuSgsLISmpiZXx8LCQu6kTwghhChbrTlnnJSUhJs3byIpKQklJSW4efMmbt68iTdv3gAAevXqBQcHB4waNQrx8fE4c+YMZsyYAX9/f+jq6gIAfHx8IBQK4efnh7/++gsHDx7E4sWLaSY1IYQQXtWaS5uCg4Oxbds2brlt27YAgHPnzqFbt25QV1fHsWPHEBAQgM6dO0MsFsPHxwcRERHcNnp6eoiJicHEiRPh6uoKfX19BAUFISgoqNr3hxBCCClTa5Lx1q1bK7zGuIyVlRWio6M/WqdVq1a4cOGCEiMjRD69evVC7qsU6ItL+Q6FEFLD1JpkTIiq7FkzGTmSdlBTU+1Zm0mTJtEjFAkh5aJkTOo9Y0NdiHRN+A6DEFKP1ZoJXIQQQkhdRcmYEEII4RkNU5N679djl5HFkqGlpYWBAweqrJ+hQ4fiRUY6TAy0Eb1llsr6IYTUPpSMSb23Ycdp7hGKqkzGeXl5yM17izwtTZX1QQipnWiYmhBCCOEZJWNCCCGEZ5SMCSGEEJ5RMiaEEEJ4RsmYEEII4RklY0IIIYRnlIwJIYQQnlEyJoQQQnhGN/0g9Z51IyNo6RjAwMBApf18++23KHnzFNpq2SrthxBS+1AyJvVe1NIJeK3bVeX9dOnShR6hSAgpFw1TE0IIITyjZEwIIYTwjIapCakm9+/fB8tNhoS9QMvmjfkOhxBSg1AyJvXerPCdeJG3D3p6eli4cKHK+pk5cyYyMjJgaqSHszvnq6wfQkjtQ8mY1HvXEv7lHqFICCF8oHPGhBBCCM8oGRNCCCE8o2RMCCGE8KzWJOOwsDB06tQJWlpaaNiwocz6W7duYfjw4bC0tIRYLIa9vT1WrVolUy8hIQHu7u4Qi8Vo1KgRFi5cCMZYNewBIYQQUr5aM4GrsLAQgwcPhpubGzZv3iyz/vr16zA2Nsb27dthaWmJP//8E+PHj4e6ujoCAwMBADk5OejZsyc8PDwQFxeHhw8fws/PDxKJBNOnT6/uXSKEEEIA1KJkHBoaCgDYunVruevHjh0rtWxra4vY2FgcOHCAS8Y7duxAfn4+tm7dCqFQCEdHRzx8+BDLly9HUFAQBAKBSveBEEIIKU+tGaaujOzsbKmb/8fGxsLd3R1CoZAr8/LyQkpKCh4/flxhOwUFBcjJyZF6EUIIIcpSZ5NxbGwsfv31V3z99ddcWVpaGkxNTaXqlS2npaVV2FZ4eDj09PS4l6WlpWqCJoQQUi/xOkwdEhLCDT9XJC4uDq6urgq1e+fOHQwYMADBwcHo2bOn1LoPh6LLJm99bIh6zpw5CAoK4pZzcnIoIdchX/Zpj8wiQ2hra6u0n927d0OjIAmiwqcq7YcQUvvwmowDAwMxbNiwj9axsbFRqM27d+/C09MT/v7+mD9f+paDZmZmMkfA6enpACBzxPw+oVAoNbRN6paAUb2q5RGKEokEmg20IGwgUnlfhJDahddkbGRkBCMjI6W1d+fOHXh6esLX1xdhYWEy693c3DB37lwUFhZCU1MTAHDq1ClYWFgolPTLjqbf5OUrJW7Cv6j9Udi3b98n63322Wcyf1vz5s3D33///cltv/zyS4wY0BFFBfR3Q0hdUJYDlHJ5LKslnjx5wuLj41loaCjT1tZm8fHxLD4+nr1+/Zoxxthff/3FjI2N2YgRI1hqair3Sk9P59p49eoVMzU1ZcOHD2cJCQnswIEDTFdXl0VERCgUS3JyMgNAL3rRi170ohd79OhRlXOcgLHacccLPz8/bNu2Tab83Llz6NatW4Xnn62traVmSickJGDixIm4evUq9PX1MWHCBAQHByt0WVNpaSlSUlKgo6Ojssuhys5LJycnQ1dXVyV9VCfan5qrLu0LQPtTk9WlfQHeXbFjZWWFrKyscm9GpYhak4zrm5ycHOjp6SE7O7tO/NHS/tRcdWlfANqfmqwu7Qug3P2ps5c2EUIIIbUFJWNCCCGEZ5SMayihUIgFCxbUmUuqaH9qrrq0LwDtT01Wl/YFUO7+0DljQgghhGd0ZEwIIYTwjJIxIYQQwjNKxoQQQgjPKBkTQgghPKNkXMsUFBTAyckJAoEAN2/e5DschT1+/Bjjxo1DkyZNIBaL0bRpUyxYsACFhYV8hya3devWoUmTJhCJRHBxccHFixf5DqlSwsPD0a5dO+jo6MDExAQDBw7EgwcP+A5LKcLDwyEQCDB16lS+Q6m0Z8+eYeTIkTA0NISWlhacnJxw/fp1vsOqlOLiYsyfP5/73Nva2mLhwoUoLS3lOzS5XLhwAd7e3rCwsIBAIMChQ4ek1jPGEBISAgsLC4jFYnTr1g137txRqA9KxrXMrFmzYGFhwXcYlXb//n2UlpZi48aNuHPnDlasWIENGzZg7ty5fIcmlz179mDq1KmYN28e4uPj0bVrV/Tp0wdJSUl8h6aw33//HRMnTsTly5cRExOD4uJi9OrVC7m5uXyHViVxcXGIjIxE69at+Q6l0rKystC5c2doaGjg+PHjuHv3LpYtW1blWy7y5ccff8SGDRuwdu1a3Lt3D0uWLMHSpUuxZs0avkOTS25uLtq0aYO1a9eWu37JkiVYvnw51q5di7i4OJiZmaFnz554/fq1/J1U+e7WpNr89ttvzM7Ojt25c4cBYPHx8XyHpBRLlixhTZo04TsMubRv355NmDBBqszOzo7Nnj2bp4iUJz09nQFgv//+O9+hVNrr16/ZZ599xmJiYpi7uzubMmUK3yFVyrfffsu6dOnCdxhK069fPzZ27Fipsi+++IKNHDmSp4gqDwA7ePAgt1xaWsrMzMzYDz/8wJXl5+czPT09tmHDBrnbpSPjWuL58+fw9/fHL7/8Ai0tLb7DUars7GwYGBjwHcYnFRYW4vr16+jVq5dUea9evfDnn3/yFJXyZGdnA0Ct+F1UZOLEiejXrx969OjBdyhVcuTIEbi6umLw4MEwMTFB27Zt8dNPP/EdVqV16dIFZ86cwcOHDwEAt27dwqVLl9C3b1+eI6u6xMREpKWlSX0vCIVCuLu7K/S9wOvzjIl8GGPw8/PDhAkT4OrqKvUUqtru0aNHWLNmDZYtW8Z3KJ/04sULlJSUwNTUVKrc1NQUaWlpPEWlHIwxBAUFoUuXLnB0dOQ7nErZvXs3bty4gbi4OL5DqbJ///0X69evR1BQEObOnYurV69i8uTJEAqFGD16NN/hKezbb79FdnY27OzsoK6ujpKSEoSFhWH48OF8h1ZlZZ/98r4Xnjx5Inc7dGTMo5CQEAgEgo++rl27hjVr1iAnJwdz5szhO+QKybsv70tJSUHv3r0xePBgfPXVVzxFrrgPH5vJGFPZozSrS2BgIG7fvo1du3bxHUqlJCcnY8qUKdi+fTtEIhHf4VRZaWkpnJ2dsXjxYrRt2xZff/01/P39sX79er5Dq5Q9e/Zg+/bt2LlzJ27cuIFt27YhIiKi3Mfi1lZV/V6gI2MeBQYGYtiwYR+tY2Njg0WLFuHy5csy9z91dXXFiBEjasQftLz7UiYlJQUeHh5wc3NDZGSkiqNTDiMjI6irq8scBaenp8v8V1ybTJo0CUeOHMGFCxfQuHFjvsOplOvXryM9PR0uLi5cWUlJCS5cuIC1a9eioKAA6urqPEaoGHNzczg4OEiV2dvbY//+/TxFVDUzZ87E7Nmzue+IVq1a4cmTJwgPD4evry/P0VWNmZkZgHdHyObm5ly5ot8LlIx5ZGRkBCMjo0/WW716NRYtWsQtp6SkwMvLC3v27EGHDh1UGaLc5N0X4N0lGx4eHnBxcUFUVBTU1GrHAI2mpiZcXFwQExODQYMGceUxMTEYMGAAj5FVDmMMkyZNwsGDB3H+/Hk0adKE75AqrXv37khISJAqGzNmDOzs7PDtt9/WqkQMAJ07d5a5zOzhw4ewtrbmKaKqycvLk/mcq6ur15pLmz6mSZMmMDMzQ0xMDNq2bQvg3fyS33//HT/++KPc7VAyrgWsrKyklrW1tQEATZs2rXVHMikpKejWrRusrKwQERGBjIwMbl3Zf5g1WVBQEEaNGgVXV1fuqD4pKQkTJkzgOzSFTZw4ETt37sThw4eho6PDHfHr6elBLBbzHJ1idHR0ZM51SyQSGBoa1spz4NOmTUOnTp2wePFiDBkyBFevXkVkZGStGUX6kLe3N8LCwmBlZYWWLVsiPj4ey5cvx9ixY/kOTS5v3rzBP//8wy0nJibi5s2bMDAwgJWVFaZOnYrFixfjs88+w2effYbFixdDS0sLPj4+8neipNnepBolJibW2kuboqKiGIByX7XFf//7X2Ztbc00NTWZs7Nzrb0UqKLfQ1RUFN+hKUVtvrSJMcaOHj3KHB0dmVAoZHZ2diwyMpLvkCotJyeHTZkyhVlZWTGRSMRsbW3ZvHnzWEFBAd+hyeXcuXPlflZ8fX0ZY+8ub1qwYAEzMzNjQqGQff755ywhIUGhPugRioQQQgjPasfJOkIIIaQOo2RMCCGE8IySMSGEEMIzSsaEEEIIzygZE0IIITyjZEwIIYTwjJIxIYQQwjNKxoQQQgjPKBkTUgsJBAIcOnSI7zCUzs/Pj3vKV1X37/0nia1cuVIp8RGiKpSMCakh3k9EGhoaMDU1Rc+ePbFlyxaZG+qnpqaiT58+crVb2xJ37969Fdq/isyYMQOpqam17v7tpH6iZExIDVKWiB4/fozjx4/Dw8MDU6ZMQf/+/VFcXMzVMzMzk3mkZl0hFAqVsn/a2towMzOrdU9sIvUTJWNCapCyRNSoUSM4Oztj7ty5OHz4MI4fP46tW7dy9d4/2i0sLERgYCDMzc0hEolgY2OD8PBwAP97hvSgQYMgEAi45UePHmHAgAEwNTWFtrY22rVrh9OnT0vFYmNjg8WLF2Ps2LHQ0dGBlZWVzFODnj59imHDhsHAwAASiQSurq64cuUKt/7o0aNwcXGBSCSCra0tQkNDpf6pkMfjx48hEAjw66+/omvXrhCLxWjXrh0ePnyIuLg4uLq6QltbG71795Z6ChghtQklY0JqOE9PT7Rp0wYHDhwod/3q1atx5MgR/Prrr3jw4AG2b9/OJd24uDgAQFRUFFJTU7nlN2/eoG/fvjh9+jTi4+Ph5eUFb29vJCUlSbW9bNkyuLq6Ij4+HgEBAfjmm29w//59rg13d3ekpKTgyJEjuHXrFmbNmsUNqZ88eRIjR47E5MmTcffuXWzcuBFbt25FWFhYpd6HBQsWYP78+bhx4wYaNGiA4cOHY9asWVi1ahUuXryIR48eITg4uFJtE8I7pT9rihBSKb6+vmzAgAHlrhs6dCizt7fnlgGwgwcPMsYYmzRpEvP09GSlpaXlbvt+3Y9xcHBga9as4Zatra3ZyJEjueXS0lJmYmLC1q9fzxhjbOPGjUxHR4dlZmaW217Xrl3Z4sWLpcp++eUXZm5uXmEM5b0HZY8M3bRpE1e2a9cuBoCdOXOGKwsPD2ctWrSQadPa2pqtWLGiwj4JqQka8PuvACFEHowxCASCctf5+fmhZ8+eaNGiBXr37o3+/fujV69eH20vNzcXoaGhiI6ORkpKCoqLi/H27VuZI+PWrVtzPwsEApiZmSE9PR0AcPPmTbRt2xYGBgbl9nH9+nXExcVJHQmXlJQgPz8feXl50NLSkmvfy4vF1NQUANCqVSupsrLYCKltKBkTUgvcu3cPTZo0KXeds7MzEhMTcfz4cZw+fRpDhgxBjx49sG/fvgrbmzlzJk6ePImIiAg0a9YMYrEYX375JQoLC6XqaWhoSC0LBAJuGFosFn805tLSUoSGhuKLL76QWScSiT66bXnej6XsH5MPyz6cdU5IbUHJmJAa7uzZs0hISMC0adMqrKOrq4uhQ4di6NCh+PLLL9G7d2+8fPkSBgYG0NDQQElJiVT9ixcvws/PD4MGDQLw7vzv48ePFYqrdevW2LRpE9fPh5ydnfHgwQM0a9ZMoXYJqY8oGRNSgxQUFCAtLQ0lJSV4/vw5Tpw4gfDwcPTv3x+jR48ud5sVK1bA3NwcTk5OUFNTw969e2FmZoaGDRsCeDcr+syZM+jcuTOEQiH09fXRrFkzHDhwAN7e3hAIBPjuu+8UPqocPnw4Fi9ejIEDByI8PBzm5uaIj4+HhYUF3NzcEBwcjP79+8PS0hKDBw+Gmpoabt++jYSEBCxatKiqbxUhdQrNpiakBjlx4gTMzc1hY2OD3r1749y5c1i9ejUOHz5c4fWy2tra+PHHH+Hq6op27drh8ePH+O2336Cm9u7jvWzZMsTExMDS0hJt27YF8C6B6+vro1OnTvD29oaXlxecnZ0VilVTUxOnTp2CiYkJ+vbti1atWuGHH37g4vTy8kJ0dDRiYmLQrl07dOzYEcuXL4e1tXUV3iFC6iYBY4zxHQQhhADvJqO9evVKqXcMs7GxwdSpUzF16lSltUmIstGRMSGkRomOjoa2tjaio6Or1M7ixYuhra0tM0OckJqIjowJITVGeno6cnJyAADm5uaQSCSVbuvly5d4+fIlAMDY2Bh6enpKiZEQVaBkTAghhPCMhqkJIYQQnlEyJoQQQnhGyZgQQgjhGSVjQgghhGeUjAkhhBCeUTImhBBCeEbJmBBCCOEZJWNCCCGEZ/8PzxNdCKAiXzoAAAAASUVORK5CYII=", "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((-5, 0), width=15, height=10, fc=\"b\", zorder=0, alpha=0.1)\n", "ax.add_patch(sky)\n", "\n", "# Aquifer:\n", "ground = plt.Rectangle(\n", " (-5, -122),\n", " width=15,\n", " height=98,\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", " (-0.5, -(122)), width=1, height=122, fc=np.array([200, 200, 200]) / 255, zorder=1\n", ")\n", "ax.add_patch(well)\n", "\n", "# Confining Unit\n", "conf = plt.Rectangle(\n", " (-5, -24),\n", " width=15,\n", " height=24,\n", " fc=np.array([100, 100, 100]) / 255,\n", " zorder=0,\n", " alpha=0.9,\n", ")\n", "ax.add_patch(conf)\n", "\n", "# Wellhead\n", "wellhead = plt.Rectangle(\n", " (-0.6, 0), width=1.2, height=4, fc=np.array([200, 200, 200]) / 255, zorder=2, ec=\"k\"\n", ")\n", "ax.add_patch(wellhead)\n", "\n", "# Screen for the well:\n", "screen = plt.Rectangle(\n", " (-0.5, -(122)),\n", " width=1,\n", " height=98,\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 = 1,y = 1.5, dx = 0, dy = 1, color = \"#00035b\")\n", "# ax.add_patch(pumping_arrow)\n", "ax.text(x=1, y=2.5, s=r\"$ Q = 10.16$ L\", fontsize=\"large\")\n", "\n", "\n", "# last line\n", "line = plt.Line2D(xdata=[-200, 1200], ydata=[0, 0], color=\"k\")\n", "ax.add_line(line)\n", "\n", "# Water table\n", "# wt = plt.Line2D(xdata= [-200,1200], ydata = [0,0], color = \"b\")\n", "# ax.add_line(wt)\n", "\n", "ax.text(0.6, -35, s=\"Ln-2\", fontsize=\"large\")\n", "# ax.text(6.9, -0.5, \"Ln-3\", fontsize = 'large')\n", "ax.set_xlim([-5, 10])\n", "ax.set_ylim([-122, 10])\n", "ax.set_xlabel(\"Distance [m]\")\n", "ax.set_ylabel(\"Relative height [m]\")\n", "ax.set_title(\"Conceptual Model - Dawsonville Example\");" ] }, { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "### Set basic parameters" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "b = 98 # aquifer thickness, m\n", "zt = -24 # top of aquifer, m\n", "zb = zt - b # bottom of aquifer, m\n", "rw = 0.076 # well radius of Ln-2 Well, m\n", "rc = 0.076 # casing radius of Ln-2 Well, m\n", "Q = 0.01016 # slug volume in m^3, m" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load data\n", "\n", "Data for the Dawsonville test is available in a text file, where the first column is the time data, in days and in the second column is the head displacement in meters" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "data = np.loadtxt(\"data/dawsonville_slug.txt\")\n", "to = data[:, 0]\n", "ho = data[:, 1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create First Model - single layer\n", "\n", "We begin with a single layer model built in ModelMaq.\n", "Details on setting up the model can be seen in: [Confined 1 - Oude Korendijk](confined1_oude_korendijk.ipynb).\n", "\n", "The slug well is set accordingly. Details on setting up the ```Well``` object can be seen in: [Slug 1 - Pratt County](slug1_pratt_county.ipynb)." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "self.neq 1\n", "solution complete\n" ] } ], "source": [ "ml = ttm.ModelMaq(kaq=10, z=[zt, zb], Saq=1e-4, tmin=1e-6, tmax=1e-3, topboundary=\"conf\")\n", "w = ttm.Well(ml, xw=0, yw=0, rw=rw, rc=rc, tsandQ=[(0, -Q)], layers=0, wbstype=\"slug\")\n", "ml.solve()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Model calibration both simultaneous wells\n", "\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": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "................................\n", "Fit succeeded.\n" ] } ], "source": [ "# unknown parameters: kay, Saq\n", "ca = ttm.Calibrate(ml)\n", "ca.set_parameter(name=\"kaq0\", initial=10, pmin=0)\n", "ca.set_parameter(name=\"Saq0\", initial=1e-4)\n", "ca.seriesinwell(name=\"obs\", element=w, t=to, h=ho)\n", "ca.fit()" ] }, { "cell_type": "code", "execution_count": 7, "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
kaq00.4209700.0183814.3662470.0inf10.0000[0.42096992918521536]
Saq00.0000170.00000531.133619-infinf0.0001[1.698434412425383e-05]
\n", "
" ], "text/plain": [ " optimal std perc_std pmin pmax initial \\\n", "kaq0 0.420970 0.018381 4.366247 0.0 inf 10.0000 \n", "Saq0 0.000017 0.000005 31.133619 -inf inf 0.0001 \n", "\n", " parray \n", "kaq0 [0.42096992918521536] \n", "Saq0 [1.698434412425383e-05] " ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "rmse: 0.004409603758298554\n" ] } ], "source": [ "display(ca.parameters)\n", "print(\"rmse:\", ca.rmse())" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdsAAAE/CAYAAAAOr2mgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVDElEQVR4nO3deVxU1fvA8c8wDPsuCKgIuO9LqInmQia5p2bx1VTcStNKv2YulevP1Mq1Rc19K9fU3DIp9/VrprllLomkoggqqAgMM/f3x8gUAsrgDMPyvF+vG8yZc+88w3F65tx77jkqRVEUhBBCCGExNtYOQAghhCjqJNkKIYQQFibJVgghhLAwSbZCCCGEhUmyFUIIISxMkq0QQghhYZJshRBCCAuTZCuEEEJYmCRbIYQQwsIk2RZyS5YsQaVSoVKp2L17d5bnFUWhQoUKqFQqmjdvbtbXVqlUjBs3zuT9oqOjUalULFmyJFf1MjYbGxs8PT1p0aIFO3bsyFvQZta8efNMf9fk5GTGjRuXbVvkl4SEBEaNGkW1atVwdnbG3d2dKlWq0KNHD06ePGmsl/FvJzo62qLx9OrVi6CgILMec9y4cahUKrMeszix9GdXZGVr7QCEebi6urJw4cIsCXXPnj1cunQJV1dX6wRmBu+++y7dunVDp9Nx7tw5xo8fT5s2bdi5cydNmza1dniZJCcnM378eACzf7nJjfv379OwYUPu37/PBx98QO3atXn48CHnz59n/fr1nDhxglq1agHQtm1bDh06hL+/f77HKURxI8m2iIiIiODbb7/l66+/xs3NzVi+cOFCQkNDSUpKsmJ0z6Zs2bI0bNgQgMaNG1OxYkWaNWvGwoULC1yytba1a9dy8eJFdu7cSVhYWKbnhg4dil6vNz728fHBx8cnv0MsFpKTk3FycrJ2GKIAkdPIRUTXrl0BWLlypbEsMTGR77//nj59+mS7z+3btxk4cCClS5fGzs6OcuXK8dFHH5GampqpXlJSEm+++SYlSpTAxcWFVq1acf78+WyPeeHCBbp160bJkiWxt7enatWqfP3112Z6lwb16tUD4ObNm5nKb9y4Qf/+/SlTpgx2dnYEBwczfvx40tPTM9WbM2cOtWvXxsXFBVdXV6pUqcKHH35ofD6nU5RPO+0aHR1tTF7jx483nv7u1asXALdu3eKtt94iICAAe3t7fHx8aNy4MT///HNe/xRZJCQkAOTYW7Wx+ecjn937ad68OTVq1ODo0aM0adIEJycnypUrx5QpUzIlaoAzZ84QHh6Ok5MTPj4+DBo0iK1bt+Z4SePfFEVh9uzZ1KlTB0dHRzw9PenSpQt//fVX3t44sHr1asLDw/H398fR0ZGqVasycuRIHjx4YKyzfPlyVCoVhw4dyrL/hAkT0Gg0XL9+3Vj2888/06JFC9zc3HBycqJx48b88ssvmfbL+Pfy22+/0aVLFzw9PSlfvnyOcWb83Xfu3Gn8XLm5udGzZ08ePHjAjRs3eP311/Hw8MDf359hw4ah1WozHaOwfnaLM0m2RYSbmxtdunRh0aJFxrKVK1diY2NDRERElvopKSmEhYWxbNkyhg4dytatW+nevTufffYZnTt3NtZTFIWOHTuyfPly3n//fTZs2EDDhg1p3bp1lmOePXuW+vXrc/r0aaZNm8aWLVto27Yt7733nvHUqjlcvnwZgEqVKhnLbty4QYMGDfjpp58YM2YMP/74I3379mXy5Mm8+eabxnqrVq1i4MCBNGvWjA0bNrBx40b++9//Zvofcl75+/uzfft2APr27cuhQ4c4dOgQo0ePBqBHjx5s3LiRMWPGsGPHDhYsWMBLL71kTJDmEBoaCkDPnj3ZuHFjno5948YN3njjDbp3786mTZto3bo1o0aNYsWKFcY6sbGxNGvWjD///JM5c+awbNky7t27xzvvvJOr1+jfvz9DhgzhpZdeYuPGjcyePZszZ87QqFGjLF+icuvChQu0adOGhQsXsn37doYMGcKaNWto3769sU5ERAR+fn5Zkkh6ejrffPMNnTp1olSpUgCsWLGC8PBw3NzcWLp0KWvWrMHLy4uXX345S8IF6Ny5MxUqVGDt2rXMnTv3qfH269cPd3d3Vq1axccff8x3333Hm2++Sdu2balduzbr1q0jMjKSadOm8eWXXxr3K8yf3WJNEYXa4sWLFUA5evSosmvXLgVQTp8+rSiKotSvX1/p1auXoiiKUr16daVZs2bG/ebOnasAypo1azId79NPP1UAZceOHYqiKMqPP/6oAMqsWbMy1fvkk08UQBk7dqyx7OWXX1bKlCmjJCYmZqr7zjvvKA4ODsrt27cVRVGUy5cvK4CyePHiJ763jHqffvqpotVqlZSUFOXEiRNKaGio4u/vr1y+fNlYt3///oqLi4ty5cqVTMeYOnWqAihnzpwxxuLh4fHE1x07dqyS3Ucj42/979dt1qxZpr/rrVu3svxdMri4uChDhgx54mubw4QJExQ7OzsFUAAlODhYGTBggPL7779nqpfT+wGUI0eOZKpbrVo15eWXXzY+/uCDDxSVSmX8u2Z4+eWXFUDZtWuXsSwyMlIJDAw0Pj506JACKNOmTcu0799//604Ojoqw4cPf+p7zKmNMuj1ekWr1Sp79uxRgEzvfezYsYqdnZ1y8+ZNY9nq1asVQNmzZ4+iKIry4MEDxcvLS2nfvn2m4+p0OqV27dpKgwYNssQyZsyYp8atKP/83d99991M5R07dlQAZfr06ZnK69Spozz33HPGx4Xhsyuykp5tEdKsWTPKly/PokWLOHXqFEePHs3xFPLOnTtxdnamS5cumcozTnlmfHPftWsXAG+88Uamet26dcv0OCUlhV9++YVOnTrh5OREenq6cWvTpg0pKSkcPnw4T+9rxIgRaDQaHBwcqFOnDqdPn2bz5s2ZRrhu2bKFsLAwSpUqlem1M77F79mzB4AGDRpw9+5dunbtyg8//EB8fHyeYsqLBg0asGTJEiZOnMjhw4eznBrMyb/fT3p6OspTlqAePXo0MTExLFq0iP79++Pi4sLcuXMJCQnJdJkhJ35+fjRo0CBTWa1atbhy5Yrx8Z49e6hRowbVqlXLVC/jcsaTbNmyBZVKRffu3TO9Lz8/P2rXrm08Ba0oSpb3/iR//fUX3bp1w8/PD7VajUajoVmzZgD88ccfxnpvv/02APPnzzeWffXVV9SsWdM4BuDgwYPcvn2byMjITK+v1+tp1aoVR48ezXI25NVXX33qe/+3du3aZXpctWpVwDBw7fHyf//tC9NnV/xDkm0RolKp6N27NytWrGDu3LlUqlSJJk2aZFs3ISEBPz+/LNcmS5Ysia2trfH0Y0JCAra2tpQoUSJTPT8/vyzHS09P58svv0Sj0WTa2rRpA5DnxDZ48GCOHj3K/v37mTp1KlqtlldeeSXTKdKbN2+yefPmLK9dvXr1TK/do0cPFi1axJUrV3j11VcpWbIkzz//PFFRUXmKzRSrV68mMjKSBQsWEBoaipeXFz179uTGjRs57hMdHZ3lPWV8cXgSX19fevfuzdy5czl58iR79uzBzs6OwYMHP3Xfx9sawN7enocPHxofJyQk4Ovrm+3rPs3NmzdRFAVfX98s7+3w4cPGtlq6dGmW53Ny//59mjRpwpEjR5g4cSK7d+/m6NGjrF+/HiBT7L6+vkRERPDNN9+g0+k4efIk+/bty3QKPONUdpcuXbLE8Omnn6IoCrdv384Ug6mjur28vDI9trOzy7E8JSXF+LgwfXbFP2Q0chHTq1cvxowZw9y5c/nkk09yrFeiRAmOHDmCoiiZPrRxcXGkp6fj7e1trJeenk5CQkKmD+3jCcLT0xO1Wk2PHj0YNGhQtq8ZHBycp/dUpkwZ46Coxo0b4+fnR/fu3Rk7dixfffUVAN7e3tSqVSvH95xxHQ6gd+/e9O7dmwcPHrB3717Gjh1Lu3btOH/+PIGBgTg4OACQmpqKvb29cb9n/R+Ot7c3M2fOZObMmcTExLBp0yZGjhxJXFyc8VpvdnEfPXo0U1nlypVNfu2mTZsSHh7Oxo0biYuLo2TJknl6DxlKlCiR7bXVJ31xyODt7Y1KpWLfvn2Z/r4ZMsrat2+f5b3nZOfOnVy/fp3du3cbe7MAd+/ezbb+4MGDWb58OT/88APbt2/Hw8MjUw8w49//l19+aRwJ/7jHv1jk132/hemzK/4hybaIKV26NB988AHnzp0jMjIyx3otWrRgzZo1bNy4kU6dOhnLly1bZnweICwsjM8++4xvv/2W9957z1jvu+++y3Q8JycnwsLCOH78OLVq1TJ+S7eEN954gwULFjB//nw++OADAgMDadeuHdu2baN8+fJ4enrm6jjOzs60bt2atLQ0OnbsyJkzZwgMDDSenj558iT169c31t+8efNTj5mRKP7dk8pO2bJleeedd/jll184cOBAjvXs7OyMXzRy4+bNm/j4+GQadQyg0+m4cOECTk5OeHh45Pp4OWnWrBlTp07l7NmzmU4lr1q16qn7tmvXjilTpnDt2jVef/31HOuVKFEi2152djKSzuPJ+5tvvsm2fkhICI0aNeLTTz/l9OnTvPXWWzg7Oxufb9y4MR4eHpw9ezbXg77yS2H+7BZnkmyLoClTpjy1Ts+ePfn666+JjIwkOjqamjVrsn//fiZNmkSbNm146aWXAAgPD6dp06YMHz6cBw8eUK9ePQ4cOMDy5cuzHHPWrFm88MILNGnShLfffpugoCDu3bvHxYsX2bx5Mzt37jTbe/z00095/vnn+b//+z8WLFjAhAkTiIqKolGjRrz33ntUrlyZlJQUoqOj2bZtG3PnzqVMmTK8+eabODo60rhxY/z9/blx4waTJ0/G3d3dmFjbtGmDl5cXffv2ZcKECdja2rJkyRL+/vvvp8bl6upKYGAgP/zwAy1atMDLywtvb288PT0JCwujW7duVKlSBVdXV44ePcr27dszjSB9VsuXL+ebb76hW7du1K9fH3d3d65evcqCBQs4c+YMY8aMMcv/TIcMGcKiRYto3bo1EyZMwNfXl++++45z584BZEn2/9a4cWPeeustevfuza+//krTpk1xdnYmNjaW/fv3U7NmTeN11dxq1KgRnp6eDBgwgLFjx6LRaPj222/5/fffc9xn8ODBREREoFKpGDhwYKbnXFxc+PLLL4mMjOT27dt06dKFkiVLcuvWLX7//Xdu3brFnDlzTIrRXAr7Z7fYsuboLPHs/j0a+UkeH42sKIqSkJCgDBgwQPH391dsbW2VwMBAZdSoUUpKSkqmenfv3lX69OmjeHh4KE5OTkrLli2Vc+fOZTvq9vLly0qfPn2U0qVLKxqNRvHx8VEaNWqkTJw4MVMdTBiN/Pnnn2f7/GuvvabY2toqFy9eVBTFMBL4vffeU4KDgxWNRqN4eXkpISEhykcffaTcv39fURRFWbp0qRIWFqb4+voqdnZ2SqlSpZTXX39dOXnyZKZj/+9//1MaNWqkODs7K6VLl1bGjh2rLFiw4KmjkRVFUX7++Welbt26ir29vQIokZGRSkpKijJgwAClVq1aipubm+Lo6KhUrlxZGTt2rPLgwYMn/h1McfbsWeX9999X6tWrp/j4+Ci2traKp6en0qxZM2X58uWZ6uY0Grl69epZjvv4iGJFUZTTp08rL730kuLg4KB4eXkpffv2VZYuXZpl9G92+yqKoixatEh5/vnnFWdnZ8XR0VEpX7680rNnT+XXX3996vvMbjTywYMHldDQUMXJyUnx8fFR+vXrp/z22285/ltLTU1V7O3tlVatWuX4Onv27FHatm2reHl5KRqNRildurTStm1bZe3atVliuXXr1lPjVpScP7M5HScyMlJxdnbOVFbQP7siK5WiPGVooxBC5NJbb73FypUrSUhIKPCnIzdv3kyHDh3YunWrcSCQEJYip5GFEHkyYcIESpUqRbly5bh//z5btmxhwYIFfPzxxwU60Z49e5YrV67w/vvvU6dOnWwneRDC3CTZCiHyRKPR8Pnnn3P16lXS09OpWLEi06dPz9XtRdY0cOBADhw4wHPPPcfSpUtl9SCRL+Q0shBCCGFhMqmFEEIIYWGSbIUQQggLk2QrhBBCWFixGyCl1+u5fv06rq6uMjBCCCGKMUVRuHfvHqVKlXriRCzmUOyS7fXr1wkICLB2GEIIIQqIv//+mzJlylj0NYpdsnV1dQUMf1w3N7cc62m1Wnbs2EF4ePgTVxsRRYu0e/Ek7V483b59m+DgYGNesKRil2wzTh27ubk9Ndk6OTnh5uYmH75iRNq9eJJ2L54y1pTOj0uKMkBKCCGEsDBJtkIIIYSFSbIVQgghLKzYXbMVQojCTK/Xk5aWZu0wCgWNRoNarbZ2GIAkWyGEKDTS0tK4fPkyer3e2qEUGh4eHvj5+Vl9XgVJtnkUm/iQy/EPCPZ2xt/d0drhCCGKOEVRiI2NRa1WExAQYPFJGAo7RVFITk4mLi4OAH9/f6vGI8k2D1YfjWHU+lPoFbBRweTONYmoX9baYQkhirD09HSSk5MpVaoUTk5O1g6nUHB0NHSE4uLiKFmypFVPKctXIxPFJj40JloAvQIfrj9NbOJD6wYmhCjSdDodAHZ2dlaOpHDJ+GKScU+ttUiyNdHl+AfGRJtBpyhExydbJyAhRLFi7WuPhU1B+XtJsjVRsLczNo+1nVqlIshbTusIIYTInlyzNZG/uyOTO1Wn0daXSFbseYg9/j7e+G5bDnbOjzaXf37au4Cd66Ofjx7bu4G9q+GnrT0UkG9eQgiR33bv3k1YWBh37tzBw8PD2uFYjCTbPIio4wPbbkFGjky4BAl5PJiNxpB4HdweJWE3cHDPfnP0ePS7Bzh6Gh5rnCRZCyFEASfJNi/U9tBvJ6TdB20ypD0w/J724J/fU+//U5Z671HZPUN5apLhMYBeCw9vG7Y8xWL3KPF6gqMXOHkZfnfyAqcSWTdnb0NClwQthBD5RpJtXqhtoUzIsx1Dr3+UgJMMSTglyfB7SuI/P1MSDeUpiZByFx7ezfxTnw66NLh/07DlOn47cPI2JF5nH8Pm4gPOJcHFF1wyfvoaErfczydEkZLf8wSkpqbywQcfsGrVKpKSkqhXrx4zZsygfv36xjoHDhzgww8/5M8//6R27dosWLCAmjVrAnDlyhXeeecd9u/fT1paGkFBQXz++ee0adPG4rGbiyRba7GxMZw6dsh5mb8nUhRDr/nhnUfbbcPP5Ee95OSMLQGS4w0/HySA9oEhQd+7btieGqfGkHRd/R5t/oafbqUebaUNP+2c8/Y+hBD5yhrzBAwfPpzvv/+epUuXEhgYyGeffcbLL7/MxYsXjXU++OADZs2ahZ+fHx9++CEdOnTg/PnzaDQaBg0aRFpaGnv37sXZ2ZmzZ8/i4uJi0ZjNTZJtYaVSPRps5QIeAbnfT/sQHsQbEvCDeLgfBw9uPfoZZ/h5P87QU35423CaO+mqYXsSB3dwDzAkX/cy4F4a3MsaYnMPMCRom4IxR6kQxVVO8wQ0reRjsR7ugwcPmDNnDkuWLKF169YAzJ8/n6ioKBYuXGjs3Y4dO5aWLVsCsHTpUsqUKcOGDRt4/fXXiYmJ4dVXXzX2dMuVK2eRWC1Jkm1xo3E0JMDcJOj0NEMCvncD7sX+8zMp1tArToqFpGuG0+EZp71vns7+WDaaR69bFjyDHtuCDYO9hBAW9aR5AiyVbC9duoRWq6Vx48bGMo1GQ4MGDfjjjz+MyTY0NNT4vJeXF5UrV+aPP/4A4L333uPtt99mx44dvPTSS7z66qvUqlXLIvFaiiRbkTNbu0e91DJPrpeSZEi6idcMPeDEa5D4N9z9GxJjDI/1Wrj9l2HLjqMXeJUzbCUqQIny//y0dzX/exOiGMqYJ+DfCdfS8wQoiuHFHp9cQlGUp044kfF8v379ePnll9m6dSs7duxg8uTJTJs2jXfffdcyQVuAJFvx7DKuPZesmv3zunRDj/juFbgbA3ei/9luXzb0nh/ehmu34dqvWfd3LQXeFcG7EvhUfrRVMQzsklHVQuSav7sjkzvX5MP1p9EpCmqVikmda1h0kFSFChWws7Nj//79dOvWDTBMnfjrr78yZMgQY73Dhw9Ttqzh2vGdO3c4f/48VapUMT4fEBDAgAEDGDBgAKNGjWL+/PmSbIXIRG375FPXqffhzmVDrzfhEty+9Oje5YuG68kZg7ku78m8n6MnlKxmSPIlq0LJ6uBbzXD9WAiRrYj6ZWlayYfo+GSCvJ0sPhrZ2dmZt99+mw8++AAvLy/Kli3LZ599RnJyMn379uX3338HYMKECZQoUQJfX18++ugjvL296dixIwBDhgyhdevWVKpUiTt37rBz506qVs3hy30BJclWWJ+9C/jVNGyPe3gH4i9C/HmI/xNunYdb5wy94od34MoBw/Zv7mXBtzr413p03FqGa8XSCxYCMPRw83Np0ClTpqDX6+nRowf37t2jXr16/PTTT3h6emaqM3jwYC5cuEDt2rXZtGmTcdEFnU7HoEGDuHr1Km5ubrRq1YoZM2bkW/zmIMlWFGyOnhBQ37D9m/YhxF+AuD8g7qxhu3n20TXjGMN2/sd/6jt4gH9tKFUHStWFUs9JAhYinzg4OPDFF1/wxRdfZHmuefPmxuu67dq1y3b/L7/80qLx5QdJtqJw0jgaeq7+j41IfHgHbp6BG6cebSch7pxhEpDLezKfinbyhtLPQel6UKYelA4BW7lfWAhhfpJsRdHi6AlBLxi2DOlphp5v7Am4fgKuHzck5OR4uLDDsD1i612JOoo/qt/vQvALhtHR0vsVQjwjSbai6LO1e3T6uA5kzLKpTTEk3Gu/wtVf4epRuHMZVfx5AjkPWx71gJ19ILARlG1k+OlbQ6avFEKYTJKtKJ40Dob5rcuEwPP9DWX3b5F+5RCX96yivF08NrEnDKOhz/5g2MBw7TfoBQhuath8qkjPVwjxVJJshcjg4oNSqTVnLyoEtWmDjUoP1357NOL5IPx9xHDt99wWwwaGeaPLNYdyYVA+zDAtpRBCPEaSrRA5sbWHwFDDBqDTGq75Ru+Fy/sg5pBhDumTqw0bgG9NqPAiVHgJyoaCWmO18IUQBYckWyFyS6355zakJu8brvte/R/8tRsu7TQk4punDNuBWYZ1g8uHQaVWUDHcsKShEKJYsvpIj9mzZxMcHIyDgwMhISHs27cvx7q7d+9GpVJl2c6dO5ePEQvxiMbBcN22xRh4azd8cAleXQi1uxoGVqUmGa71bnwbPq8AC182JOH4i089tBCiaLFqz3b16tUMGTKE2bNn07hxY7755htat27N2bNnjXNkZufPP//Eze2fdWB9fHzyI1whnsy5BNTsYtj0esMtRhd+gj9/NNzv+/dhwxY1BnyqQrUOULWDYbYrGWQlRJFm1Z7t9OnT6du3L/369aNq1arMnDmTgIAA5syZ88T9SpYsiZ+fn3FTq2WdVFHA2NgYRjqHfQgD9sF/z0CbqYaBVDa2cOsP2PMpzG0MX9WDnRMNM2AJUUw1b94808IERY3VerZpaWkcO3aMkSNHZioPDw/n4MGDT9y3bt26pKSkUK1aNT7++GPCwsJyrJuamkpqaqrxcVJSEmBYdUKr1ea4X8ZzT6ojih6LtbuTL9TtZdge3kV14Sdszm1G9dcuVAkXYe/nsPdzFO/K6Kt3Rl+jC3gEmjcGkaPC8HnXarUoioJer0ev11s7nFx78cUXqV27tnEu4927d9OiRQsSEhLw8PAw1lu3bh0ajcbs702v16MoClqtNkvHLD/b22rJNj4+Hp1Oh6+vb6ZyX19fbty4ke0+/v7+zJs3j5CQEFJTU1m+fDktWrRg9+7dNG3aNNt9Jk+ezPjx47OU79ixAyenp6/hGBUVlYt3I4oay7e7K7h0w7Z6J3wTj1P67hFKJp1CHf8n6j2TUe+ZTIJzRa56NuKaZ0O0Mo1kvijIn3dbW1v8/Py4f/8+aWlp1g4n19LT00lLSzN2dJKTkwG4d+8eNv+aIMbW1hZFUYz1zCUtLY2HDx+yd+9e0tPTMz2XEUt+UCkZM0Dns+vXr1O6dGkOHjxIaGiosfyTTz5h+fLluR701L59e1QqFZs2bcr2+ex6tgEBAcTHx2e67vs4rVZLVFQULVu2RKOR2zeKC6u2e0oiqj+3YXNmHarLe1Fh+GgqajuUSq3R14xAKf+i4TS0MKvC8HlPSUnh77//JigoCAcHB2uHkyu9e/dm2bJlT6zTs2dPFi9enKUHXK5cOfr27cv58+fZsGEDJUqUYObMmTRq1Ig333yTnTt3EhwczMKFC6lXr16Ox09JSSE6OpqAgIAsf7eEhAT8/f1JTEx8Yj4wB6t9ar29vVGr1Vl6sXFxcVl6u0/SsGFDVqxYkePz9vb22NvbZynXaDS5+lDltp4oWqzS7hpvqNfTsCXFwunv4feVqG6eRvXHD9j88QO4+EHdN6BuD/AKzt/4ioGC/HnX6XSoVCpsbGwMPUJFAW3+9cwy0TjlalDfF198wYULF6hRowYTJkxAp9Nx+PBhunTpYhzo6ujoaOzhZry/DDNnzmTSpEmMGTOGGTNmEBkZSePGjenTpw9Tp05lxIgR9OrVizNnzqDKIR4bGxtUKlW2bZufbW21ZGtnZ0dISAhRUVF06tTJWB4VFcUrr7yS6+McP34cf39/S4QohPW4+UOjdwxb7En4faVh4oz7N2DfNMMW3Azq9YEqbWXyjOJImwyTSlnntT+8DnZPv7Th7u6OnZ0dTk5O+PkZZlcrUaIEYBjo+u9rttlp06YN/fsbplMdM2YMc+bMoX79+rz22msAjBgxgtDQUG7evGk8fkFl1fNRQ4cOpUePHtSrV4/Q0FDmzZtHTEwMAwYMAGDUqFFcu3bNeBpi5syZBAUFUb16ddLS0lixYgXff/8933//vTXfhhCWlbGU4Evj4c9t8NsywyQaGUsGuvjBcz0hpBe4l7Z2tEKYTa1a/yyhmXHGs2bNmlnK4uLiJNk+SUREBAkJCUyYMIHY2Fhq1KjBtm3bCAw0jMKMjY0lJibGWD8tLY1hw4Zx7do1HB0dqV69Olu3bqVNmzbWegtC5B9bO6je0bDdjYFjSw2J9/4N2PuZobdbrQM8PwACnpd7d4s6jZOhh2mt186Pl/nXad6M08TZlRWG0dlWH2kxcOBABg4cmO1zS5YsyfR4+PDhDB8+PB+iEqKA8ygLLUZDsxHw51b43wK4sh/ObDBs/rWh4SCo0VlOMRdVKlWuTuVam52dHTqdLtNjIFNZcWD16RqFEM/A1g6qd4LeW2HAAcPpZFsHiP0dNrwFM2vB/pnw8K61IxXFVFBQEEeOHCE6Opr4+HgCAwNRqVRs2bKFW7ducf/+fWuHmC8k2QpRVPjVgA5fwtA/4MXRhuX/7l2Hn8fCjBqGaSLvZX8PuxCWMmzYMNRqNdWqVcPHxwetVsv48eMZOXIkvr6+vPPOO9YOMV9Y/TSyEMLMnLyg6TBo9C6cWgsHvzJMD3lgFhyeC3W6wQtDwDPI2pGKYqBSpUocOnQoU9no0aMZPXp0prLdu3dnehwdHZ3lWI9PCxEUFJSlrKCSnq0QRZWtPdTtDm8fhK6roEwD0KXCscXwxXOwcRAkXLJ2lEIUC5JshSjqbGygcmvouwN6bYPyLUDRwYkVhkUQ1veXpCuEhUmyFaK4UKkgqDH0WA/9foGKL4Oih5Or4Kv6sOldwy1FQgizk2QrRHFUph68scaw6H3FcENP97dlhtPLP46AB/HWjlCIIkWSrRDFWam68MZa6BtlmP5Rr4Ujc2FWbdj9KaQWj9syCpPCMiCooCgofy9JtkIICGgAkZugx0bwrwNp92H3JHSz6nBx+1fE3pGka20Za7EWpuX1CoKMZfSsvcCE3PojhPhH+TBDD/fsRu5tG4trcgwVDn/Enwfn8mfDj2ne5j/WjrDYsrW1xcnJiVu3bqHRaDKtjiOyUhSF5ORk4uLi8PDwyLJwfH6TZCuEyMzGhtiA1jS/o6abTRSDbddT2eZvKv+vPym31uLQ7jMoUd7aURY7KpUKf39/Ll++zJUrV6wdTqHh4eFRIBYpkGQrhMjicvwDUhVbFutas17XhHdtNxCp3oHD5Z/h6+chdJBh4gx7V2uHWqzY2dlRsWJFOZWcSxqNxuo92gySbIUQWQR7O2OjAr0CibgwMb0Hq/QvsaXSFhyid8GBmfD7Kmg1Cap3lhWG8pGNjQ0ODg7WDkOYKFfJdtOmTSYfuGXLljg6Opq8nxDC+vzdHZncuSYfrj+NTlFQq1S82Skch3p94fxPsH0k3LkM6/oYbhlqMw28K1g7bCEKrFwl244dO5p0UJVKxYULFyhXrlxeYhJCFAAR9cvStJIP0fHJBHk74e/+6Mtz5VZQrrlhruV90+Cv3TAnFF74LzR53zBNpBAik1wPZ7tx4wZ6vT5Xm5NT/iwsLISwLH93R0LLl/gn0WbQOEDzETDwkGH6R10a7PkU5r4AVw5lfzAhirFcJdvIyEiTTgl3794dNze3PAclhCgkSpSH7t9Dl8XgXBLiz8PiVrB5CKQkWTs6IQqMXCXbxYsX4+qa+1GHc+bMwdvbO89BCSEKEZUKanSGd/5nWLweDCsLzQ6Fiz9bNzYhCgi5K1oIYR6OnobF6yO3GNbKTboKK16FH96BlERrRyeEVZl8609KSgpffvklu3btIi4uDr1en+n53377zWzBCSEKoeAmhjV0f/k/wzzLx5fDpZ3QcbZhYJUQxZDJybZPnz5ERUXRpUsXGjRogErurxNCPM7OGVpPgWodYONAw21Cy16B59+Gl8aCRm4LFMWLycl269atbNu2jcaNG1siHiFEURLYCAbsh6jR8OsiODLH0Mt9dT7417Z2dELkG5Ov2ZYuXdqkwVJCiGLO3gXazYBua8HFF+L/hPkt4OCX8NhlKCGKKpOT7bRp0xgxYoRMhC2EME2lcHj7EFRpZ1g3d8fHsKIz3Lth7ciEsDiTk229evVISUmhXLlyuLq64uXllWkTQogcOZeAiBWGnq6tI/y1C+Y0hou/WDsyISzK5Gu2Xbt25dq1a0yaNAlfX18ZICWEMI1KBfX6QNlG8H1fuHnacItQk/eh+ShQy/ooougx+V/1wYMHOXToELVry+AGIcQzKFkF+v0M20cZJsHYNxWuHIQui8DN39rRCWFWJp9GrlKlCg8fPrRELEKI4kbjCO1nwqsLwc4FYg7CN03g8j5rRyaEWZmcbKdMmcL777/P7t27SUhIICkpKdMmhBAmq9kF+u8F3xrw4JbhntwDs0BRrB2ZEGZh8mnkVq1aAdCiRYtM5YqioFKp0Ol05olMCFG8lCgPfaNgy3/h5CqIGgN//w86zQV7ud1QFG4mJ9tdu3ZZIg4hhAA7J0NyDagPP46Ec1tgwUvwn+8MyViIQsrkZNusWTOzBjB79mw+//xzYmNjqV69OjNnzqRJkyZP3e/AgQM0a9aMGjVqcOLECbPGJISwIpUK6vcD/zqwujvcOgfzwqDLQqjY0trRCZEnubpme/LkySwLDjzJmTNnSE9Pf2q91atXM2TIED766COOHz9OkyZNaN26NTExMU/cLzExkZ49e2Y5lS2EKELK1IO3dkOZBpCaCN++BvtnynVcUSjlKtnWrVuXhISEXB80NDT0qQkTYPr06fTt25d+/fpRtWpVZs6cSUBAAHPmzHnifv3796dbt26EhobmOiYhRCHk6ge9tsBzkYACP481LGyQnmrtyIQwSa5OIyuKwujRo3FycsrVQdPS0nJV59ixY4wcOTJTeXh4OAcPHsxxv8WLF3Pp0iVWrFjBxIkTcxWPEKIQs7WH9rMMI5W3j4Dfv4Pbf8F/vgVnb2tHJ0Su5CrZNm3alD///DPXBw0NDcXR8clLaMXHx6PT6fD19c1U7uvry40b2c+VeuHCBUaOHMm+ffuwtc3d5ebU1FRSU//5Fpxxe5JWq0Wr1ea4X8ZzT6ojih5p9wLsud6oPIJQr++L6u/DKPPCSI9YCT6Vn/nQ0u7FU362d64y1u7duy0WwOPTPWbcQvQ4nU5Ht27dGD9+PJUqVcr18SdPnsz48eOzlO/YsSNXPfWoqKhcv5YoOqTdCy59mQ9pfmUGnokxsLAl/wt+j3jXamY5trR78ZKcnJxvr6VSFOuMNkhLS8PJyYm1a9fSqVMnY/ngwYM5ceIEe/bsyVT/7t27eHp6olarjWV6vR5FUVCr1ezYsYMXX3wxy+tk17MNCAggPj4eNze3HOPTarVERUXRsmVLNBrNs7xVUYhIuxdsa49d5eMfzuKm3GOB3TTq2ZxHsbFF13YmSq3/5Pm40u7FU0JCAv7+/iQmJj4xH5iD1Wb8trOzIyQkhKioqEzJNioqildeeSVLfTc3N06dOpWpbPbs2ezcuZN169YRHByc7evY29tjb2+fpVyj0eTqQ5XbeqJokXYveGITH/LxD2fRK3AXV95I+5Cpmm9ozyFsN78DSVeh+UjDrUN5JO1evORnW1t1eY2hQ4fSo0cP6tWrR2hoKPPmzSMmJoYBAwYAMGrUKK5du8ayZcuwsbGhRo0amfYvWbIkDg4OWcqFEEXP5fgH6P91Hi4VO97TDqJu7TqUOT0H9kyBe9eh7QxZOUgUOFb9FxkREUFCQgITJkwgNjaWGjVqsG3bNgIDAwGIjY3N1S1EQoiiL9jbGRsVmRKujUqNuuVYCKwI24bBb8vgfhx0WWyYjUqIAsLkhQj27t2b7YQV6enp7N271+QABg4cSHR0NKmpqRw7doymTZsan1uyZMkTB2eNGzdOZo8Sopjwd3dkcueaqB+dJlarVEzqXAN/d0eo3xdeXw62DnB+OyxtD8m3rRyxEP8wuWcbFhZGbGwsJUuWzFSemJhIWFiYLEQghLCYiPplaVrJh+j4ZIK8nQyJNkPVdtDzB/guAq79CotaQY8N4F7aegEL8YjJPducbs1JSEjA2dnZLEEJIURO/N0dCS1fInOizVC2IfT5CVxLQfyfsOhliL+Y/0EK8Zhc92w7d+4MGO6L7dWrV6YRvjqdjpMnT9KoUSPzRyiEEKYoWQX6/gTLOsLtS4aE22M9+Ne2dmSiGMt1z9bd3R13d3cURcHV1dX42N3dHT8/P9566y1WrFhhyViFECJ3PMoaerh+tSA5Hpa0g5gj1o5KFGO57tkuXrwYgKCgIIYNGyanjIUQBZuLj2ERg+/+AzEHYXkn6LoSypl3mVAhcsPka7Zjx46VRCuEKBwc3KH791D+RdA+MCzTd/4na0cliiGTk+3Nmzfp0aMHpUqVwtbWFrVanWkTQogCxc4Juq6Cym1Blwqr3oCzm6wdlShmTL71p1evXsTExDB69Gj8/f2zHZkshBAFiq09vL4UNgyA0+tgbS94dQHU6GztyEQxYXKy3b9/P/v27aNOnToWCEcIISxErYHO8ww/f18J3/cFvQ5qvWbtyEQxYHKyDQgIwEoLBQkhxLOxUcMrXxt+Hl8BG94CRQfVXrV2ZKKIM/ma7cyZMxk5ciTR0dEWCEcIISzMRg3tv4TnIkHRw4YBqE6tsXZUoogzuWcbERFBcnIy5cuXx8nJKcsSRbdvy3ykQogCzsYG2s0ElQ0cW4x68zuULtsfaGPtyEQRZXKynTlzpgXCEEKIfGZjA22ng6JD9dsyQq7MRXc2BGrLNVxhfiYn28jISEvEIYQQ+c/GBtrNQp+ejs3J71Bv7A8ae6jWwdqRiSLG5Gu2AJcuXeLjjz+ma9euxMXFAbB9+3bOnDlj1uCEEMLibGzQtZ1BjFdjVIoO/bo+3D6x2dpRiSLG5GS7Z88eatasyZEjR1i/fj33798H4OTJk4wdO9bsAQohhMXZqJlj/xabdQ2x0Wtx3tCbXT/KoClhPiYn25EjRzJx4kSioqKws7MzloeFhXHo0CGzBieEEPkhNjGFVZdt+a92IDt0IdirtDx/+B3iz+6xdmiiiDA52Z46dYpOnTplKffx8SEhIcEsQQkhRH66kpCMgop0bHlH+x57dLVwUqXisaEbXD9u7fBEEWBysvXw8CA2NjZL+fHjxyldurRZghJCiPwUWMIJFYbJetLQ0F/7Xw7rq2KrvQ8rXoVb560coSjsTE623bp1Y8SIEdy4cQOVSoVer+fAgQMMGzaMnj17WiJGIYSwKH93ByLK6bF5NNW7VuXAtVaLoFRdSE6A5R3hboxVYxSFm8m3/nzyySf06tWL0qVLoygK1apVQ6fT0a1bNz7++GNLxCiEEBYX6qswsHNTriWmEeTthL+7I9T+Hha3hvg/YdkrhgXpXUpaO1RRCJmcbDUaDd9++y0TJkzg+PHj6PV66tatS8WKFS0RnxBC5Bt/dwfKerv+U+BcAnpuhEUvw+2/YHln6L3VsE6uECYwOdlmKF++POXLlzdnLEIIUfC4lYIeG2FRK7h5ClZ2MyxIr3GwdmSiEDE52SqKwrp169i1axdxcXHo9fpMz69fv95swQkhRIFQojx0XwdL2sGV/Ybl+V5bCuo891dEMWPyAKnBgwfTo0cPLl++jIuLC+7u7pk2IYQokvxrw3++A7U9nNsCW/8LstyoyCWTv5atWLGC9evX06aNrI4hhChmgptAl4Wwpif8tgxc/SHsQ2tHJQoBk3u27u7ulCtXzhKxCCFEwVe1vWG1IIA9n8Kvi6wbjygUTE6248aNY/z48Tx8+NAS8QghRMFXrzc0G2H4fev7cG6bdeMRBZ7Jp5Ffe+01Vq5cScmSJQkKCsqyePxvv/1mtuCEEKLAaj4Kkq7D8eWwrg9EboKABtaOShRQJifbXr16cezYMbp3746vry8qlcoScQkhRMGmUkG7mXD/JlzYAd9FQL+fDSOXhXiMycl269at/PTTT7zwwguWiEcIIQoPtS28tgSWtIXrx0lf1pnfwtcQUKasYQYqIR4x+ZptQEAAbm5ulohFCCEKHztn6LaG+46lsU2MxnZ1N16csp3VR2UuZfEPk5PttGnTGD58ONHR0WYJYPbs2QQHB+Pg4EBISAj79u3Lse7+/ftp3LgxJUqUwNHRkSpVqjBjxgyzxCGEEHkVq3OlU+IQ7irOPGdzkRm2X/Px+pPEJspAUmFg8mnk7t27k5ycTPny5XFycsoyQOr27du5Ptbq1asZMmQIs2fPpnHjxnzzzTe0bt2as2fPUrZs2Sz1nZ2deeedd6hVqxbOzs7s37+f/v374+zszFtvvWXqWxFCCLO4HP+AC/rSvJU2lOV2k2mlPspw5Vui40PldLIA8pBsZ86cabYXnz59On379qVfv37GY//000/MmTOHyZMnZ6lft25d6tata3wcFBTE+vXr2bdvnyRbIYTVBHs7Y6OC/ylV+UA7gC/svuJN220kXlsL5QdYOzxRAJicbCMjI83ywmlpaRw7doyRI0dmKg8PD+fgwYO5Osbx48c5ePAgEydOzLFOamoqqampxsdJSUkAaLVatFptjvtlPPekOqLokXYvnp613b2dbJn4SjU+/uEsm/SNCEq/wVDbdbjt+pD0ksEo5V80Z7jCTPLzc56nWbQvXbrE4sWLuXTpErNmzaJkyZJs376dgIAAqlevnqtjxMfHo9Pp8PX1zVTu6+vLjRs3nrhvmTJluHXrFunp6YwbN87YM87O5MmTGT9+fJbyHTt24OTk9NQ4o6KinlpHFD3S7sXTs7S7MzC2LtxKUVHCvj0xN2Mpe/sAypqe7Ks0mnuOAeYLVJhFcnJyvr2Wycl2z549tG7dmsaNG7N3714++eQTSpYsycmTJ1mwYAHr1q0z6XiP36erKMpT793dt28f9+/f5/Dhw4wcOZIKFSrQtWvXbOuOGjWKoUOHGh8nJSUREBBAeHj4E0dVa7VaoqKiaNmyZZbr0qLoknYvnizS7ukvoV/ZBU3MIcKuzyG99w5ZeL6ASUhIyLfXMjnZjhw5kokTJzJ06FBcXf9ZZDksLIxZs2bl+jje3t6o1eosvdi4uLgsvd3HBQcHA1CzZk1u3rzJuHHjcky29vb22NvbZynXaDS5+lDltp4oWqTdiyeztrtGY1glaMFLqG5fQvN9JERukXVwC5D8/IybfOvPqVOn6NSpU5ZyHx8fk74l2NnZERISkuW0TVRUFI0aNcr1cRRFyXRNVgghCgwnL+i2Bhw84OpR+GGQLMtXTJncs/Xw8CA2NtbYu8xw/PhxSpcubdKxhg4dSo8ePahXrx6hoaHMmzePmJgYBgwwjN4bNWoU165dY9myZQB8/fXXlC1blipVqgCG+26nTp3Ku+++a+rbEEKI/OFdAV5fBis6w+l14FMZmg23dlQin5mcbLt168aIESNYu3YtKpUKvV7PgQMHGDZsGD179jTpWBERESQkJDBhwgRiY2OpUaMG27ZtIzAwEIDY2FhiYv6ZhUWv1zNq1CguX76Mra0t5cuXZ8qUKfTv39/UtyGEEPmnXDNoOw02D4Zdn0CJClCjM7GJD7kc/4Bgb2e5H7eIUymKaec0tFotvXr1YtWqVSiKgq2tLTqdjm7durFkyRLUarWlYjWLpKQk3N3dSUxMfOoAqW3bttGmTRu5dleMSLsXT/nW7j99BIe+AltHdjRczIBf9OgVsFHB5M41iaifdTIfYTkJCQl4e3s/NR+Yg8nXbDUaDd9++y0XLlxgzZo1rFixgnPnzrF8+fICn2iFEMKqWk6ACi0h/SE19w2khHIHAL0CH64/LdM7FmF5us8WoFy5cpQrV86csQghRNFmo4YuC0me8yL+iReZbzeNiLQxpGKHTlGIjk+W08lFlMk92y5dujBlypQs5Z9//jmvvfaaWYISQogiy8Gd+51XcEdxoY7NX3yqmQcoqFUqgryfPtGOKJxMTrZ79uyhbdu2WcpbtWrF3r17zRKUEEIUZSUDq3Ki4Sy0ipqO6oO8bbuFSZ1rSK+2CDM52d6/fx87O7ss5RqNxjjvsBBCiCcLa92F5BaTABhuu4oIt7NWjkhYksnJtkaNGqxevTpL+apVq6hWrZpZghJCiOLAvekAqNcHFQp83w/izlk7JGEhJg+QGj16NK+++iqXLl3ixRcNK1n88ssvrFy5krVr15o9QCGEKNJafQq3zsOV/bCqK/T7xTDzlChSTO7ZdujQgY0bN3Lx4kUGDhzI+++/z9WrV/n555/p2LGjBUIUQogizNYOXl8K7mXh9l+wrg/o0q0dlTCzPN3607Zt22wHSQkhhMgDZ2/ouhIWtoS/dsHPY+HlT6wdlTAjk3u2QgghLMCvBnScbfj90Fdwco114xFmZXKy1el0TJ06lQYNGuDn54eXl1emTQghRB5V7wRNhhl+3/QuXD8OQGziQw5eipcZpgoxk5Pt+PHjmT59Oq+//jqJiYkMHTqUzp07Y2Njw7hx4ywQohBCFCNhH0GlVpCeAqveYOP+32g8ZSfd5h+h8ZSdrD4a8/RjiALH5GT77bffMn/+fIYNG4atrS1du3ZlwYIFjBkzhsOHD1siRiGEKD5sbKDzPChREZKuUWrH26gVw4ApmUO58DI52d64cYOaNWsC4OLiQmJiIgDt2rVj69at5o1OCCGKIwd36LqSdI0LDWzO8bHtcuNTGXMoi8LF5GRbpkwZYmNjAahQoQI7duwA4OjRo9jb25s3OiGEKK68K5LU2jBgKtI2itfUuwFkDuVCyuRk26lTJ3755RcABg8ezOjRo6lYsSI9e/akT58+Zg9QCCGKK6/nXuFUpUEATLRdxHM2l2QO5ULK5Pts/73iT5cuXShTpgwHDx6kQoUKdOjQwazBCSFEcVfzPxNJ+fYvHC79yFrPr1FXibB2SCIP8ryebYaGDRvSsGFDc8QihBDicTY2OLw+H+a3QB3/J6yNhJ6bDDNPiUIjV8l206ZNuT6g9G6FEMLM7F3hP9/B/DCIOQQ/fQhtp1o7KmGCXCXb3M55rFKp0Ol0zxKPEEKI7HhXgM7zYWUEHJ0PpepA3e7WjkrkUq4GSOn1+lxtkmiFEMKCKreC5h8aft/yX7h2zLrxiFyTuZGFEKIwafoBVG4LujRY3QPu37J2RCIX8pRsf/nlF9q1a0f58uWpUKEC7dq14+effzZ3bEIIIR5nYwOd5hpnmGJtL2JvJ8ncyQWcycn2q6++olWrVri6ujJ48GDee+893NzcaNOmDV999ZUlYhRCCPFvDm6GAVN2rnBlP9unvylzJxdwJifbyZMnM2PGDFauXMl7773He++9x3fffceMGTOYNGmSJWIUQgjxOJ9K3H75CwB6226no81+mTu5ADM52SYlJdGqVass5eHh4SQlJZklKCGEEE93zqMps9I7ATBFM5/qqmiZO7mAMjnZdujQgQ0bNmQp/+GHH2jfvr1ZghJCCPF0wd7OfKF7lZ26OjiotHxjNx1v1T2ZO7kAMnkGqapVq/LJJ5+we/duQkNDATh8+DAHDhzg/fff54svvjDWfe+998wXqRBCiEz83R2Z1Lk2768fxAbVxwTZ3GRLqUX4uXSxdmjiMSYn24ULF+Lp6cnZs2c5e/assdzDw4OFCxcaH6tUKkm2QghhYRH1y9K0kg83L5RB/9Or+CUcgV/GQ/j/WTs08S8mJ9vLly9bIg4hhBB55O/uiH+9RuDwNazrDQe/gFJ1oUZna4cmHnnmSS10Oh0nTpzgzp075ohHCCFEXtXoDI0HG37/YRDcPANAbOJDuQ/XykxOtkOGDDGeLtbpdDRt2pTnnnuOgIAAdu/ebXIAs2fPJjg4GAcHB0JCQti3b1+OddevX0/Lli3x8fHBzc2N0NBQfvrpJ5NfUwghiqwWY6FcGGiTYVU31h88ReMpO+U+XCszOdmuW7eO2rVrA7B582aio6M5d+4cQ4YM4aOPPjLpWKtXrzbud/z4cZo0aULr1q2Jicn+H8PevXtp2bIl27Zt49ixY4SFhdG+fXuOHz9u6tsQQoiiyUYNXRaBR1m4E43Xj4NA0QPIfbhWZHKyjY+Px8/PD4Bt27bx2muvUalSJfr27cupU6dMOtb06dPp27cv/fr1o2rVqsycOZOAgADmzJmTbf2ZM2cyfPhw6tevT8WKFZk0aRIVK1Zk8+bNpr4NIYQoupy8IOJbdGoHmqt/Z6jtWuNTch+udZg8QMrX15ezZ8/i7+/P9u3bmT17NgDJycmo1epcHyctLY1jx44xcuTITOXh4eEcPHgwV8fQ6/Xcu3cPLy+vHOukpqaSmppqfJwx8YZWq0Wr1ea4X8ZzT6ojih5p9+KpSLa7d1USX/wc76h3ecf2B07pg/lJ3wAbFZR2tyta7zWP8vNvYHKy7d27N6+//jr+/v6oVCpatmwJwJEjR6hSpUqujxMfH49Op8PX1zdTua+vLzdu3MjVMaZNm8aDBw94/fXXc6wzefJkxo8fn6V8x44dODk9/cbvqKioXMUiihZp9+Kp6LW7O17OrWny4Eema+bQOc2PusGlOX5gJ3LxzdBJzC8mJ9tx48ZRo0YN/v77b1577TXs7e0BUKvVWXqpuaFSqTI9VhQlS1l2Vq5cybhx4/jhhx8oWbJkjvVGjRrF0KFDjY+TkpIICAggPDwcNze3HPfTarVERUXRsmVLNBpNLt6JKAqk3YunIt3urcJJWd4F56v72VJyDvznF3D0sHZUBUJCQkK+vZbJyRagS5ess5NERkaadAxvb2/UanWWXmxcXFyW3u7jVq9eTd++fVm7di0vvfTSE+va29sbvxD8m0ajydWHKrf1RNEi7V48Fc1216DpugzmNUeTeAU2DYBuawwDqYq5/GzrXCXbL774grfeegsHB4dM0zFmJ7ezRtnZ2RESEkJUVBSdOnUylkdFRfHKK6/kuN/KlSvp06cPK1eupG3btrl6LSGEKNacS8B/VsDCl+Hiz7Dz/+ClccQmPuRy/AOCvZ3xd3e0dpRFWq6S7YwZM3jjjTdwcHBgxowZOdYzdYrGoUOH0qNHD+rVq0doaCjz5s0jJiaGAQMGAIZTwNeuXWPZsmWAIdH27NmTWbNm0bBhQ2Ov2NHREXd391y/rhBCFDv+taHDl7C+H+yfwcEHpeh+uDR6BWxUMLlzTSLql7V2lEVWrpLtv6doNOd0jRERESQkJDBhwgRiY2OpUaMG27ZtIzAwEIDY2NhM99x+8803pKenM2jQIAYNGmQsj4yMZMmSJWaLSwghiqRar8GN3+Hgl9T57WMqM44/CDTef9u0ko/0cC0kT9dszWngwIEMHDgw2+ceT6B5maFKCCHEv7w0njuXT+AZu4/5dtNonzqRO7gZ77+VZGsZuUq2/x7N+zTTp0/PczBCCCEszEZNasf5RM9uRpDqJl9rvqCndiSKSiPr4FpQrpLt49MhHjt2DJ1OR+XKlQE4f/48arWakJAQ80cohBDCrPx8/fmx2Vy893SlkfosY1iBfYdp0qu1oFwl2127dhl/nz59Oq6urixduhRPT08A7ty5Q+/evWnSpIllohRCCGFWrV98kdsec2FTL3qqd4DNTqCXtcMqskyeG3natGlMnjzZmGgBPD09mThxItOmTTNrcEIIISzH67lOEPax4cHWYXDlkHUDKsJMTrZJSUncvHkzS3lcXBz37t0zS1BCCCHySdNhUK0j6LWwujvcjZH1by3A5NHInTp1onfv3kybNo2GDRsCcPjwYT744AM6d+5s9gCFEEJYkEoFHWfD7Utw4xR3FnWh5a0R3Fcc5P5bMzK5Zzt37lzatm1L9+7dCQwMJDAwkDfeeIPWrVsbVwASQghRiNg5w39WonPyxjPpT6bazkGFXta/NSOTk62TkxOzZ88mISGB48eP89tvv3H79m1mz56Ns7OzJWIUQghhaR4BnG0ym1TFllbqowy1XQfI+rfmYnKyzeDs7EytWrWoXbu2JFkhhCgCvKs15cP0fgC8a7uRDjYHUKtUcv+tGeQ52QohhCha/N0dadBxEN+ktwfgc808vgnTyf23ZiDJVgghhFFE/bJ0eH8ut8u0wF6l5aXf/wt3Y56+o3giSbZCCCEy8fd0wavHUvCtAQ9uwXf/gZQka4dVqEmyFUIIkZW9K3RbDS6+EHeGlNW9OHjhpoxMziNJtkIIIbLnXga6riTdxgGHy7/w59J3aTxlJ6uPymllU0myFUIIkaNYl2oMTu0PQG/bn4i0+VHuvc0DSbZCCCFydDn+AVt1zzNZ2xWA0bYreFF1VO69NZEkWyGEEDkK9nbGRgXf6NrxbXoLbFQKX2i+omL6nzKHsglMnhtZCCFE8eHv7sjkzjX5cP1pxqT3ooxNPM1sfof13Xk5aTQxSkmZQzkXpGcrhBDiiSLql2X/yDBWvNmYyu9+j9anOo6pCSzWfIonSTKHci5IshVCCPFU/u6OhJYvgZ+PDyeazueaUoLyNrEssJuGA6kyh/JTSLIVQghhkjJly9FbO4JExYkQmwt8qfkKO5UeJzsbuYabA7lmK4QQwiT+7o707dSatzY8YJlmEi3Vx1jivYZOs1XoFZVcw82G9GyFEEKYLKJ+WWaOeJvoZjNRUNHoziYGqw3L8sk13Kwk2QohhMgTf3dHKr/Ynb8ajAdgsO0GItU/AbIO7uMk2QohhHgmTo3fYkZ6FwDGa5bSweagrIP7GEm2Qgghnom/uyOlOoxhmS4cgGmaOSx64Y6sg/svkmyFEEI8s4gGgbR8fwnxQe3RqHQ0Oz4ULu+zdlgFhiRbIYQQZuHv4Yx3j8VQqRWkp8DK/8DVX60dVoEgyVYIIYT5qDXw2lIIbgpp99Eu68yti8esHZXVSbIVQghhXhoHvq/0Ocf0FdGkJaJa3pEfd+4q1gsXyKQWQgghzCo28SEfbLqEizKcb+0+oaZNNPX2RNL154+5qC9dLCe9sHrPdvbs2QQHB+Pg4EBISAj79uV8QT02NpZu3bpRuXJlbGxsGDJkSP4FKoQQIlcuxz9Ar0ASzvRIG8VZfSA+qkS+1XxCsCq2WE56YdVku3r1aoYMGcJHH33E8ePHadKkCa1btyYmJibb+qmpqfj4+PDRRx9Ru3btfI5WCCFEbmSsgQtwF1feSBvFH/oAfFV3WWk3kSBVbLGb9MKqyXb69On07duXfv36UbVqVWbOnElAQABz5szJtn5QUBCzZs2iZ8+euLu753O0QgghciNjDVy1ypBxE3Gje9pH/Kkvg5/qDqvt/o9KNteL1aQXVrtmm5aWxrFjxxg5cmSm8vDwcA4ePGilqIQQQphDRP2yNK3kQ3R8MkHeTuw9f4se6z9mqeYTqtr8zSaXSTg8DAX3GtYONV9YLdnGx8ej0+nw9fXNVO7r68uNGzfM9jqpqamkpqYaHyclJQGg1WrRarU57pfx3JPqiKJH2r14kna3DG8nW7zLugHQuY4/ocHtuB5bl7S9/XC4dQplaTvSu64Df+tcFszP9rb6aGTVo9MMGRRFyVL2LCZPnsz48eOzlO/YsQMnp6efwoiKijJbLKLwkHYvnqTd80eU30Aa3J+Kz8NLKEvac6T8f7ntUjnf40hOzr9rxlZLtt7e3qjV6iy92Li4uCy93WcxatQohg4danyclJREQEAA4eHhuLm55bifVqslKiqKli1botFozBaPKNik3Ysnaff8tfbYVVr8Oor5mqk8zzlCL03ldpt5/OkaSmAJJ/zdHfIljoSEhHx5HbBisrWzsyMkJISoqCg6depkLI+KiuKVV14x2+vY29tjb2+fpVyj0eTqQ5XbeqJokXYvnqTdLS828SEf/3AWveJEZNoIvtZ8QQuO47W5N+u0A9isvJBv9+DmZ1tb9TTy0KFD6dGjB/Xq1SM0NJR58+YRExPDgAEDAEOv9Nq1ayxbtsy4z4kTJwC4f/8+t27d4sSJE9jZ2VGtWjVrvAUhhBAmyLgHFyAFe/pr/8tnzKOzej+z7GZTQnuPD9eraFrJp0itGmTVZBsREUFCQgITJkwgNjaWGjVqsG3bNgIDAwHDJBaP33Nbt25d4+/Hjh3ju+++IzAwkOjo6PwMXQghRB5k3IObkXDTseV97QASFWd62/7EGM1ySqfHs/X3SrStXbrIJFyrzyA1cOBAoqOjSU1N5dixYzRt2tT43JIlS9i9e3em+oqiZNkk0QohROHw+D24No/+Oz69J1O0/wGgr+2PlIp6mxenbGf10ewnOSpsrD4aWQghRPGS3T24H64/zVxdB64r3nyumUsb9f8oqbrLwPVDi8QpZUm2Qggh8p2/u6MxgWYk360nY5m4FW6meTLPbhr1bM6z3u5jDh3wJLRxs0KdcK1+GlkIIYTwd3ekbS1/bFRwRKlK57TxXNb7UkYVz8tHejL2s88K9SllSbZCCCEKhH9fz72klKZj2v+xX1cdZ1Uq8zTTubpxAlt+v1ooVwuSZCuEEKLAiKhflv0jw/i4bVUScaGXdgRL01sC8L5mLfbretBqyuZC18uVZCuEEKJA+fcp5XRsGZvemw+0b5GqaGip/o1Nmo9YtmFroerhSrIVQghR4Dx+i9BaXXM6p43jb70PgTZxfK8ZzfX931o5ytyTZCuEEKJAyjil/FXXutio4IwSTLu0T9ilq40tOj47cLfQnE6WW3+EEEIUWP7ujrSr7ciDtHRGfX+KRFzoo/2AuqqL/KZU4tf1pwvFfbjSsxVCCFHgRdQvyxfdDNP1Ktjwm1IJAJ2iEB2ff0vl5ZUkWyGEEIVCSKAnNo8td65WqQjyfvra5NYmyVYIIUSh8PigKbVKxaTONQr8KWSQa7ZCCCEKkcfnVS4MiRYk2QohhChk/j2vcmEhp5GFEEIIC5NkK4QQQliYJFshhBDCwiTZCiGEEBZW7AZIKYoCQFJS0hPrabVakpOTSUpKQqPR5EdoogCQdi+epN2Lp3v37gH/5AVLKnbJNuOPGxAQYOVIhBBCFAQJCQm4u7tb9DVUSn6k9AJEr9dz/fp1XF1dUT26Mbp+/focPXo0U72kpCQCAgL4+++/cXNzs0aoOcou3oJwXFP3z239p9XL6/PS7uY5rrS7ZRXEds/LvvnR7qY+l5iYSNmyZblz5w4eHh5Pje1ZFLuerY2NDWXKlMlUplarc/yAubm5FbgP35PiteZxTd0/t/WfVi+vz0u7m+e40u6WVRDbPS/75ke75/U5GxvLD1+SAVLAoEGDrB2CSSwV77Me19T9c1v/afXy+ry0u3mOK+1uWQWx3fOyb360e16fyw/F7jRybiUlJeHu7k5iYmKB+6YrLEfavXiSdi+e8rPdpWebA3t7e8aOHYu9vb21QxH5SNq9eJJ2L57ys92lZyuEEEJYmPRshRBCCAuTZCuEEEJYmCRbIYQQwsIk2QohhBAWJslWCCGEsDBJtmZw+fJlwsLCqFatGjVr1uTBgwfWDknkA1tbW+rUqUOdOnXo16+ftcMR+Sg5OZnAwECGDRtm7VBEPrh37x7169enTp061KxZk/nz55t8DLn1xwyaNWvGxIkTadKkCbdv38bNzQ1b22I3E2ax4+3tTXx8vLXDEFbw0UcfceHCBcqWLcvUqVOtHY6wMJ1OR2pqKk5OTiQnJ1OjRg2OHj1KiRIlcn0M6dk+ozNnzqDRaGjSpAkAXl5ekmiFKMIuXLjAuXPnaNOmjbVDEflErVbj5OQEQEpKCjqdzuRl+Yp8st27dy/t27enVKlSqFQqNm7cmKXO7NmzCQ4OxsHBgZCQEPbt25fr41+4cAEXFxc6dOjAc889x6RJk8wYvcgrS7c7GKZ6CwkJ4YUXXmDPnj1milw8i/xo92HDhjF58mQzRSzMIT/a/e7du9SuXZsyZcowfPhwvL29Tdq/yHfBHjx4QO3atenduzevvvpqludXr17NkCFDmD17No0bN+abb76hdevWnD17lrJlywIQEhJCampqln137NiBVqtl3759nDhxgpIlS9KqVSvq169Py5YtLf7eRM4s3e6lSpUiOjqaUqVKcfr0adq2bcupU6dkXl0rs3S7Hz16lEqVKlGpUiUOHjxo8fcjcic/Pu8eHh78/vvv3Lx5k86dO9OlSxd8fX1zH6RSjADKhg0bMpU1aNBAGTBgQKayKlWqKCNHjszVMQ8ePKi8/PLLxsefffaZ8tlnnz1zrMJ8LNHuj2vVqpVy9OjRvIYoLMAS7T5y5EilTJkySmBgoFKiRAnFzc1NGT9+vLlCFmaQH5/3AQMGKGvWrDFpnyJ/GvlJ0tLSOHbsGOHh4ZnKw8PDc/2ttX79+ty8eZM7d+6g1+vZu3cvVatWtUS4wkzM0e537twxfgu+evUqZ8+epVy5cmaPVZiPOdp98uTJ/P3330RHRzN16lTefPNNxowZY4lwhZmYo91v3rxJUlISYLh8tHfvXipXrmxSHEX+NPKTxMfHo9PpspwK8PX15caNG7k6hq2tLZMmTaJp06YoikJ4eDjt2rWzRLjCTMzR7n/88Qf9+/fHxsYGlUrFrFmz8PLyskS4wkzM0e6i8DFHu1+9epW+ffuiKAqKovDOO+9Qq1Ytk+Io1sk2g0qlyvRYUZQsZU/SunVrWrdube6whIU9S7s3atSIU6dOWSIsYWHP+nnP0KtXLzNFJPLDs7R7SEgIJ06ceKbXL9ankb29vVGr1Vm+3cTFxZl24VsUKtLuxZO0e/FUUNq9WCdbOzs7QkJCiIqKylQeFRVFo0aNrBSVsDRp9+JJ2r14KijtXuRPI9+/f5+LFy8aH1++fJkTJ07g5eVF2bJlGTp0KD169KBevXqEhoYyb948YmJiGDBggBWjFs9K2r14knYvngpFu+dp3HMhsmvXLgXIskVGRhrrfP3110pgYKBiZ2enPPfcc8qePXusF7AwC2n34knavXgqDO0ucyMLIYQQFlasr9kKIYQQ+UGSrRBCCGFhkmyFEEIIC5NkK4QQQliYJFshhBDCwiTZCiGEEBYmyVYIIYSwMEm2QgghhIVJshWigNu9ezcqlYq7d+/m+2urVCpUKhUeHh5PrDdu3Djq1KmT6XHGvjNnzrRojEIUBpJshShAmjdvzpAhQzKVNWrUiNjYWNzd3a0S0+LFizl//rxJ+wwbNozY2FjKlCljoaiEKFyK/EIEQhR2dnZ2+Pn5We31PTw8KFmypEn7uLi44OLiglqttlBUQhQu0rMVooDo1asXe/bsYdasWcZTsNHR0VlOIy9ZsgQPDw+2bNlC5cqVcXJyokuXLjx48IClS5cSFBSEp6cn7777Ljqdznj8tLQ0hg8fTunSpXF2dub5559n9+7deYp1ypQp+Pr64urqSt++fUlJSTHDX0CIokuSrRAFxKxZswgNDeXNN98kNjaW2NhYAgICsq2bnJzMF198wapVq9i+fTu7d++mc+fObNu2jW3btrF8+XLmzZvHunXrjPv07t2bAwcOsGrVKk6ePMlrr71Gq1atuHDhgklxrlmzhrFjx/LJJ5/w66+/4u/vz+zZs5/pvQtR1MlpZCEKCHd3d+zs7HBycnrqaWOtVsucOXMoX748AF26dGH58uXcvHkTFxcXqlWrRlhYGLt27SIiIoJLly6xcuVKrl69SqlSpQDDddXt27ezePFiJk2alOs4Z86cSZ8+fejXrx8AEydO5Oeff5berRBPID1bIQohJycnY6IF8PX1JSgoCBcXl0xlcXFxAPz2228oikKlSpWM11NdXFzYs2cPly5dMum1//jjD0JDQzOVPf5YCJGZ9GyFKIQ0Gk2mxyqVKtsyvV4PgF6vR61Wc+zYsSyDlv6doIUQliHJVogCxM7OLtOgJnOpW7cuOp2OuLg4mjRp8kzHqlq1KocPH6Znz57GssOHDz9riEIUaZJshShAgoKCOHLkCNHR0bi4uODl5WWW41aqVIk33niDnj17Mm3aNOrWrUt8fDw7d+6kZs2atGnTJtfHGjx4MJGRkdSrV48XXniBb7/9ljNnzlCuXDmzxCpEUSTXbIUoQIYNG4ZaraZatWr4+PgQExNjtmMvXryYnj178v7771O5cmU6dOjAkSNHchzxnJOIiAjGjBnDiBEjCAkJ4cqVK7z99ttmi1OIokilKIpi7SCEEAWTSqViw4YNdOzYMU/7BwUFMWTIkCyzYglR3EjPVgjxRF27djV52sVJkybh4uJi1p65EIWZ9GyFEDm6ePEiAGq1muDg4Fzvd/v2bW7fvg2Aj4+P1eZ1FqKgkGQrhBBCWJicRhZCCCEsTJKtEEIIYWGSbIUQQggLk2QrhBBCWJgkWyGEEMLCJNkKIYQQFibJVgghhLAwSbZCCCGEhUmyFUIIISzs/wFfzlVutzlGqAAAAABJRU5ErkJggg==", "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)\n", "plt.semilogx(to, ho, \".\", label=\"obs\")\n", "plt.semilogx(tm, hm[0], label=\"ttim\")\n", "plt.xlabel(\"time [d]\")\n", "plt.ylabel(\"displacement [m]\")\n", "plt.title(\"Model Results - Single-layer model\")\n", "plt.legend()\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In general, the single-layer model seems to be performing well, with a good visual fit between observations and the model." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Analysis and comparison of simulated values\n", "\n", "We now compare the values in TTim and add the results of the modelling done in MLU by Yang (2020)." ] }, { "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
MLU0.4133000.0000190.004264
ttim0.4209700.0000170.004410
\n" ], "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ta = pd.DataFrame(\n", " columns=[\"k [m/d]\", \"Ss [1/m]\"],\n", " index=[\"MLU\", \"ttim\"],\n", ")\n", "ta.loc[\"MLU\"] = [0.4133, 1.9388e-05]\n", "ta.loc[\"ttim\"] = ca.parameters[\"optimal\"].values\n", "ta[\"RMSE\"] = [0.004264, ca.rmse()]\n", "ta.style.set_caption(\"Comparison of parameter values and error under different models\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Results are similar between both models. The RMSE of MLU is slightly better than the one from TTim." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### References\n", "\n", "* Cooper Jr, H.H., Bredehoeft, J.D., Papadopulos, I.S., 1967. Response of a finite diameter well to an instantaneous charge of water. Water Resources Research 3, 263–269\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 }