{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"# 3. Leaky Aquifer Test - Texas Hill"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Import packages"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"import timflow.transient as tft\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",
"This example concerns a pumping test conducted in the location of Texas Hill and is taken from the AQTESOLV examples (Duffield, 2007). A pumping well is screened in an aquifer located between 20 ft and 70 ft depths. The aquifer is overlain by an aquitard. The formation at the base of the aquifer is considered impermeable.\n",
"\n",
"Three observation wells are located at 40, 80, and 160 ft distance. They are called OW1, OW2, and OW3, respectively. Pumping lasted for 420 minutes at a rate of 4488 gallons per minute. "
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"### Load data"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"# data of OW1\n",
"data1 = np.loadtxt(\"data/texas40.txt\")\n",
"t1 = data1[:, 0] # days\n",
"h1 = data1[:, 1] # meters\n",
"r1 = 40 * 0.3048 # distance between obs1 to pumping well in m\n",
"\n",
"# data of OW2\n",
"data2 = np.loadtxt(\"data/texas80.txt\")\n",
"t2 = data2[:, 0]\n",
"h2 = data2[:, 1]\n",
"r2 = 80 * 0.3048 # distance between obs2 to pumping well in m\n",
"\n",
"# data of OW3\n",
"data3 = np.loadtxt(\"data/texas160.txt\")\n",
"t3 = data3[:, 0]\n",
"h3 = data3[:, 1]\n",
"r3 = 160 * 0.3048 # distance between obs3 to pumping well in m"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"### Parameters and model"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"# known parameters\n",
"Q = (4488 * 0.00378541) * 60 * 24 # discharge, m^3/d\n",
"b1 = 20 * 0.3048 # overlying aquitard thickness, m\n",
"b2 = 50 * 0.3048 # aquifer thickness, m\n",
"zt = -b1 # top of aquifer, m\n",
"zb = -b1 - b2 # bottom of aquifer, m\n",
"rw = 0.5 * 0.3048 # well radius, m"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"The overlying layer is modeld as an aquitard without storage (```Sll```, the storage parameter, is set to zero in ModelMaq)."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"self.neq 1\n",
"solution complete\n"
]
}
],
"source": [
"ml = tft.ModelMaq(\n",
" kaq=10,\n",
" z=[0, zt, zb],\n",
" Saq=0.001,\n",
" Sll=0,\n",
" c=10,\n",
" tmin=0.001,\n",
" tmax=1,\n",
" topboundary=\"semi\",\n",
")\n",
"w = tft.Well(ml, xw=0, yw=0, rw=rw, tsandQ=[(0, Q)], layers=0)\n",
"ml.solve()"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"### Estimate aquifer parameters\n",
"The hydraulic parameters of the aquifer (`kaq` and `Saq`) and the aquitard (`c`) are estimated using observations in all three observation wells. "
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"...................................................\n",
"Fit succeeded.\n"
]
}
],
"source": [
"# unknown parameters: kaq, Saq, c\n",
"cal = tft.Calibrate(ml)\n",
"cal.set_parameter(name=\"kaq\", layers=0, initial=10)\n",
"cal.set_parameter(name=\"Saq\", layers=0, initial=1e-4)\n",
"cal.set_parameter(name=\"c\", layers=0, initial=100)\n",
"cal.series(name=\"OW1\", x=r1, y=0, t=t1, h=h1, layer=0)\n",
"cal.series(name=\"OW2\", x=r2, y=0, t=t2, h=h2, layer=0)\n",
"cal.series(name=\"OW3\", x=r3, y=0, t=t3, h=h3, layer=0)\n",
"cal.fit(report=False)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" optimal | \n",
"
\n",
" \n",
" \n",
" \n",
" | kaq_0_0 | \n",
" 224.627650 | \n",
"
\n",
" \n",
" | Saq_0_0 | \n",
" 0.000213 | \n",
"
\n",
" \n",
" | c_0_0 | \n",
" 43.880381 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" optimal\n",
"kaq_0_0 224.627650\n",
"Saq_0_0 0.000213\n",
"c_0_0 43.880381"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"RMSE: 0.060 m\n"
]
}
],
"source": [
"display(cal.parameters.loc[:, [\"optimal\"]])\n",
"print(f\"RMSE: {cal.rmse():.3f} m\")"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAm8AAAEqCAYAAABDS7fKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAACJJklEQVR4nOzdeVxU5f7A8c+ZYd9BQEBBBBfAFRVlMQE3bHNJ08pyyXa7Lfdm2XbNftlyb93KvGV1LTVNzbIsM9cY9wX3BVAUERSQfd9n5vfHKEqCgjKy+H2/XvMyzpw58xw6M+fL93me76Po9Xo9QgghhBCiRVA1dQOEEEIIIUT9SfAmhBBCCNGCSPAmhBBCCNGCSPAmhBBCCNGCSPAmhBBCCNGCSPAmhBBCCNGCSPAmhBBCCNGCmDR1A24lnU5Hamoqtra2KIrS1M0RQgjRBPR6PYWFhXh4eKBSSQ5DtDy3VfCWmpqKp6dnUzdDCCFEM5CSkkL79u2buhlCNNhtFbzZ2toChg+snZ1dE7dGtGSVlZVs2LCB4cOHY2pq2tTNEa2MXF/GVVBQgKenZ/U9QYiW5rYK3i51ldrZ2UnwJm5KZWUlVlZW2NnZyc1VNDq5vm4NGT4jWirp7BdCCCGEaEEkeBNCCCGEaEEkeBNCCCGEaEFuqzFvQgghRFPQ6XRUVFQ0dTNEM2Vqaopara73/hK8CSGEEEZUUVHBmTNn0Ol0Td0U0Yw5ODjg5uZWr4k0Erw1UFFuGXkZpTi4WmLjaNHUzRFCCNGM6fV60tLSUKvVeHp6SlFgcRW9Xk9JSQkZGRkAuLu7X/c1Erw1QOyOVDRL4tHrQVEg4mE/AsI8mrpZQgghmqmqqipKSkrw8PDAysqqqZsjmilLS0sAMjIycHV1vW4XqvwJUE9FuWXVgRuAXg+apfEU5ZY1bcOEEEI0W1qtFgAzM7Mmbolo7i4F95WVldfdV4K3esrLKK0O3C7R6yA/o7RpGiSEEKLFkILA4noaco1I8FZPDq6WKAroqjLQVpxCr69CUYG9q2VTN00IIYQQtxEJ3urJxtGCiIf90FYcoLL4V8rzv8TZLYbctAT0MoNICCGEELeITFhogIAwDzJOd+P41lTKCvNIOb6NlOPbsG3jgt/AcPwHRuDi5d3UzRRCCCFEKyaZtwaKmDSZp79axLg33qFbxFDMLK0ozM4kZvWPLJ7xLItnPMve1T9SkJXZ1E0VQgghbkpKSgqPPvooHh4emJmZ0aFDB55//nmys7MBmDlzJn5+fjVeEx8fj6IoTJkypcb2hQsXYm5uTmmpYaz4nDlzCA0NxcrKCgcHh1txOq2GBG83QKVS06FHb0Y8/QJPffUd97wwE99+wajUJmQmJ7Ht+4V8/eyjrJg9kyOb11NWXNTUTRZCCNHCpeWXsvN0Fmn5t2aiXGJiIv369SMhIYFly5Zx6tQp5s+fz+bNmwkJCSEnJ4fIyEhOnDhBenp69euio6Px9PREo9HUOF50dDTBwcHVZTEqKiq4//77efrpp2/J+bQm0m16k0zNzOkaMpCuIQMpLSokYfcO4rZrOBd3jHOxhsef33xBx8Ag/O+IwCcwCBOZMi6EEKIBVsQk8+qqo+j0oFLgvft6MCHIy6jvOX36dMzMzNiwYUN1wOXl5UVgYCC+vr68/vrrfPjhh5iamqLRaHjggQcA0Gg0TJ8+nTlz5pCUlIS3t3f19qlTp1Yff/bs2YAhIycaRjJvjcjSxpaeQ0cw4a33eXzeNwx8cDJt2nuhrariVMwufvvPe8x/8hHWz59L8rEjMtFBCCHEdaXll1YHbgA6Pby26phRM3A5OTmsX7+eZ555pjpwu8TNzY2JEyeyYsUKrKysCAoKIjo6uvp5jUbDkCFDCAsLq96emJhIcnIykZGRRmvz7UQyb0Zi5+LKgNH303/UODLPniFuu4b4HVsoysnmWPQGjkVvwKaNM36hgwwTHTp0lDpAQgghrnImq7g6cLtEq9eTlFWCu71xylUlJCSg1+vx9/ev9Xl/f39yc3PJzMwkMjKSlStXAhAbG0tZWRmBgYEMGjSoOtum0WiwsLAgODjYKO293UjwZmSKouDq7YOrtw+DHprCubhjxG7TkLBnB0XZWez7bRX7fltFm/Ze+N8Rif/AcOycXZu62UIIIZqJjs7WqBRqBHBqRcHb2fjLben/Wp2+FhEREcyZM4e0tDQ0Gg0DBw5ErVYTHh7O/PnzAUM2LjQ0FHNzc2M3+bYg3aa3kKJS4dmtJ1FPPcdTX37HyL+/Ruf+oahNTMg+l8z2ZYv4evqjLJ/1Ckc2raO0qLCpmyyEEKKJudtb8t59PVBf7J1RKwrv3tfdaFk3gE6dOqEoCnFxcbU+HxcXh6OjIy4uLoSFhWFmZkZ0dDTR0dGEh4cDEBQURFZWFomJiWg0GgYPHmy09t5uJPPWREzMzOg8IJTOA0IpKyri5J4dxG/XkBJ3jPPxxzkff5zN38ynY2A//AdG4NM3CFMz+YtFCCFuRxOCvBjUxYWkrBK8na2MGrgBtGnThmHDhvH555/z4osv1hj3lp6eztKlS5k0aRKKomBpacmAAQPQaDRs2bKFGTNmAGBqakpwcDALFiwgJSVFxrs1IgnemgELGxt6Domi55AoCrIyid+xhfjtGjKTkzi9bzen9+3GzNKKzgNC8R8YgWe3HqhU6qZuthBCiFvI3d7S6EHblebNm0doaChRUVG88847dOzYkePHjzNjxgzatWvHnDlzqveNjIzk448/BqBPnz7V28PDw/nwww+xtrYmKCioxvGTk5PJyckhOTkZrVbLoUOHAEPWz8bGxvgn2IJJt2kzY+fsQv9R45j073lM/vc8+o8ah20bFypKSziu2cSP77zB189MRfPdAi6cOV2v8QhCCCFEQ3Xu3Jl9+/bh4+PD+PHj8fX15YknniAyMpJdu3bh5ORUvW9kZCSFhYWEhYVhYnI5LxQeHk5hYSEDBw7E1NS0xvH/+c9/EhgYyKxZsygqKiIwMJDAwED27dt3y86xpVL0t9Hdv6CgAHt7e/Lz87Gzs2vq5tSbXqfjfHwscds1nNy9vUbRX6d2nvgPjMB/YDj2rm51HqMot4y8jFIcXC2xcbS4Fc1u1SorK1m7di133XXXVV9IQtwsub6M61beC8rKyjhz5gwdO3bEwkK+e0XdGnKtSLdpC6CoVLQP6E77gO5ETn2SM4f2Eb9Nw+kDe8k5n8KOFd+xY8V3eHQNwH9gBF1DBmJpe/kLKXZHKpol8ej1oCgQ8bAfAWEeTXdCQgghhLhhEry1MCampnQOCqFzUAjlJcUk7NlJ3PZoko8fJfVELKknYole+CXevfviPzCCtr69qgM3AL0eNEvj8QpwkgycEEII0QJJ8NaCmVtZ0z1yGN0jh1GYk8WJHVuJ276FjKTTJO7fS+L+vZiYWaDDB7WZPyoTTxRFhV4H+RmlErwJIYQQLZAEb62ErZMz/e69j3733kf2uWTitmuI276FgswLQCy6ilhQrFGbd8fEojv2rrduxpIQQgghGo8Eb61Qm/ZeDHxgEmETHiH1RBw7fvydlGN7QF+MtmwP2rK9rJ9/iJ5DovDp0x+1iVwGQgghREshd+1WTFEU2vkFMP6NAPIzC4nduoOzR7ZyPv4ISYf2k3RoP9aOTnSPGEaPwcOxd23b1E0WQgghxHVI8HabsHexJWTsCELGjiA3PZWjf27guGYTxbk57Pl5BXt++QHvnoH0HDICn76SjRNCCCGaK7lD34Yc3TwY9NAUwsZP5PS+PRzZvJ6zRw6SdPgASYcPYO3gSPfIS9m4umvHCSGEEOLWkxUWbmNqE1O6BA9k3Ov/x7RPv6b/qHFY2TtQnJfLnp9/4H/PPc6Pc94kYc9OtFVVTd1cIYQQTUij0aAoCnl5eTd1nJKSEsaOHYudnV318by9vfnkk08apZ23AwneBAAObu7c8dAUnvj8W+79+6t06BkIej1njxzk1/+8y1fPTGHbskXkXUhv6qYKIYQwsoiICF544YUa20JDQ0lLS8Pe3v6mjr1o0SK2bdvGzp07G+V4NyolJYVHH30UDw8PzMzM6NChA88//zzZ2dnV+8ycORM/P78ar4uPj0dRFKZMmVJj+8KFCzE3N6e0tBSAOXPmEBoaipWVFQ4ODo3adgneRA1qE1O6DAgzZOPm/o/+o+/Hyt6Bkvw89v6ykgXPPcaPc97k5O7tko0TQojbiJmZGW5ubiiKclPHOX36NP7+/nTv3r1RjncjEhMT6devHwkJCSxbtoxTp04xf/58Nm/eTEhICDk5OYBhzdYTJ06Qnn45cREdHY2npycajabGMaOjowkODsbS0lCKq6Kigvvvv5+nn3660dsvwZuok0NbN+54cDJPfL6QkX9/zZCNA84eOchvH79vyMZ9v5C89LQmbqkQQojGMmXKFLZs2cKnn36KoigoikJSUtJV3aYLFy7EwcGBNWvW0LVrV6ysrBg3bhwlJSUsWrQIb29vHB0dee6559BqtYAho/fRRx+xdetWFEUhIiKi1jYkJyczatQobGxssLOzY/z48Vy4cAGA/Px81Gp19QL2Op0OJycngoODq1+/ZMkSPD096zzH6dOnY2ZmxoYNGwgPD8fLy4s777yTTZs2cf78eV5//XUABg4ciKmpaY1ATaPRMH36dHJyckhKSqqxPTIysvrn2bNn8+KLL9KjR496/+7rS4I3cV1qExM6Dwhl3Ov/x2Of/Y8BY8Zj7eBoyMat/pEFzz/OynfeuJiNq6Qot4xzJ3Ipyi1r6qYLIUTzotdDRXHTPC6tk3gdn376KSEhITz++OOkpaWRlpZWZyBUUlLC3LlzWb58OevWrUOj0TBmzBjWrl3L2rVr+e677/jyyy/58ccfAVi1ahWPP/44ISEhpKWlsWrVqquOqdPpGDVqFDk5OWzZsoWNGzeSmJjIhAkTALC3t6d3797VAdXRo0dRFIWDBw9SVFQEwJYtWwgPD6+1zTk5Oaxfv55nnnmmOkt2iZubGxMnTmTFihXo9Xqsra0JCgoiOjq6eh+NRsOQIUMICwur3p6YmEhycnKN4M2YZLapaBB7VzcGPjCJkHEPkXhgL0c2ryfp8AGSjx4i+eghzCxt0en9UJn1QG3iQMTDfgSEeTR1s4UQonmoLIF3m+g78bVUMLO+7m729vaYmZlhZWWFm9u1Kw5UVlbyxRdf4OvrC8C4ceP47rvvuHDhAjY2NgQEBBAZGUl0dDQTJkzAyckJKyur6i7Y2mzevJmjR49y5syZ6qBx8eLFdOvWjZiYGIKCgoiIiECj0fDSSy+h0WgYNmwY8fHxbN++nREjRqDRaHj55ZdrPX5CQgJ6vR5/f/9an/f39yc3N5fMzExcXV2JjIxk5cqVAMTGxlJWVkZgYCCDBg1Co9EwdepUNBoNFhYWNbJ/xiSZN3FD1CYmdO4fythXZ/PY3K8ZMGYCVnYOVJQWUlUWQ0XBN5QXrGLzt2spyC5u6uYKIYQwAisrq+rADaBt27Z4e3tjY2NTY1tGRka9jxkXF4enp2eNbF9AQAAODg7ExcUBEB4ezvbt29FqtWzZsoWIiIjqgC41NZVTp07V2SV7ib6emciIiAhOnjxJWloaGo2GgQMHolarCQ8Pr87+aTQaQkNDMTc3r/d53gzJvImbZsjGPYJXzxH88uHPaMsPo6s6i64qiYrCJJa+tp0+d95Nj8HDsbJrmllFQgjRLJhaGTJgTfXejX1IU9MaPyuKUus2nU7XqO87aNAgCgsLOXDgAFu3buXdd9/Fzc2N999/n169euHh4UHnzp1rfW2nTp1QFIW4uDjGjBlz1fNxcXE4Ojri4uICQFhYGGZmZkRHRxMdHV3dHRsUFERWVhaJiYloNBqefPLJRj3Ha5HgrYHSi9NJLkjGy84LN2spYHslJ3dbTMw7oTbrhE6bh7b8MNqK45TkZbF92SJ2rVxKl5A76D38Ltw7+zXJDCMhhGhSilKvrsumZmZmVj3J4Fbz9/cnJSWFlJSU6uxbbGwseXl5BAQEAODg4EDPnj2ZN28epqam+Pn54erqyoQJE1izZk2d490A2rRpw7Bhw/j888958cUXa4x7S09PZ+nSpUyaNKn6HmVpacmAAQPQaDRs2bKFGTNmAIbANTg4mAULFpCSknLLxruBdJs2yKqEVUT9FMW0DdOI+imKVQlXD7S8ndk4WhDxsB+KClRqB8xswhn+9IdEPf0Cbr6d0VZVEbctmmVvzmDJzBc4snk9lWUyqUEIIZobb29v9uzZQ1JSEllZWY2eObuWoUOH0qNHDyZOnMiBAwfYu3cvkyZNIjw8nH79+lXvFxERwdKlS6sDNScnJ/z9/VmxYsU1gzeAefPmUV5eTlRUFFu3biUlJYV169YxbNgw2rVrx5w5c2rsHxkZyfLlyykrK6NPnz7V28PDw/nss8+qJzZcKTk5mUOHDpGcnIxWq+XQoUMcOnSoelLFzWgxwZsxi93VR3pxOrN3zUanN1zAOr2O2btmk14sRWuvFBDmwaQ5oYx+MZBJc0LpEe5N94ihTHz3YybO+Q/dwodiYmpGRtJpNn71GV8+PZnoRV+Tk3q+qZsuhBDiopdeegm1Wk1AQAAuLi4kJyffsvdWFIXVq1fj6OjIoEGDGDp0KD4+PqxYsaLGfuHh4Wi12hpj2yIiIq7aVpvOnTuzb98+fHx8GD9+PL6+vjzxxBNERkaya9cunJycauwfGRlJYWEhYWFhmFyx9nd4eDiFhYXVJUWu9M9//pPAwEBmzZpFUVERgYGBBAYGVpc4uRmKvr4j9prYrFmzcHBw4Ny5cyxYsOCGlucoKCjA3t6e/Px87OzsGvTavWl7mbZh2lXbv4n6hiC3oFpeIepSWljAMc0mDm9cS/4VKzZ06BlIr+F34dunPyq1uglbeH2VlZWsXbuWu+6666oPrBA3S64v47qZe0FDlZWVcebMGTp27IiFhYVR30u0bA25VlrMmLfZs2cDhqKATcHLzguVoqrOvAGoFBWetnUXARS1s7S1I+je++h392iSjhzk0IbfSTwQw9kjBzl75CA2bZzpNWQEPYZEYe3g2NTNFUIIIZqVFhO83Yjy8nLKy8urfy4oKAAMf9VWVlY26FhtzNrwRv83mLP7/9CiQ6VS80b/N2hj1qbBxxKXte/Wk/bdelKQmcHRP9dzXLOJouwsdvywhF0/LadTUDA9h92Jexf/ZjXB4dL/c/l/L4xBri/jkt+raOladfD23nvvVWfsrrRhwwasrBo+ZdoMM96LH4Lj3n1kDwyhylrF2hNrG6OpAsC2DR53jaU4+Qz5CbGUZWVwcvd2Tu7ejql9Gxy6BmDr3QmVSfO5bDdu3NjUTRCtmFxfxlFSUtLUTRDipjTpXXDmzJl88MEH19wnLi4OPz+/Gzr+q6++yt///vfqnwsKCvD09GT48OE3PM4hZeEiylOzsP/hN9SaPdg/+AD299+P2l7qlzW2vb/GELP6V7Tl8VTmZ5O5dxsFxw4QEDGUHkOicGjr3mRtq6ysZOPGjQwbNkzGJIlGJ9eXcV3qhRGipWrS4O0f//gHU6ZMueY+Pj4+N3x8c3PzWqsdm5qa3vAXotc3C8hb8QO5S5dSlZFBzqdzyf3yK+xHj8Jp0mTMfTrecHvFZUW5ZRzeVIap1XBMLO5AW3EcbflhykvyObh2NQf/+BWfwH70jroH756BKKqmmTh9M9eSENcj15dxyO9UtHRNGry5uLhUVzBuKUwcHXF+6knaPDqVgj/+IHvRIspj48hbvoK85SuwCQ/HaeoUrAYMaFZjtFqavIzS6jWUFZUlJhb9UJv3YcDdapKPa0g6tJ/EAzEkHojBwc2d3sPvoVvEECysba59YCGEEKKFaz6Dh64jOTmZnJycGsXuwLDMxZVrqN0qipkZ9qNGYTdyJCV7Y8hZtIii6GiKtmyhaMsWzLt2xWnyZOzuuRuVmdktb19L5+BqiaLAlYVsVGoV/neEEDQykty08xzasJbjmk3kpaehWfw121csJuCOSHpH3YOLl3eTtV0IIYQwphZT523KlCksWrToqu3R0dHXLcZ3ibFr+1QkJZGz+Dvyfv4ZfWkpAGpnZxwfehDHBx7A5C9F/8S1xe5IRbM0Hr0OFBVETPQjIMyjxj4VZaXEbdNwaP0aslLOVm937xxAxz5D6DZoIHbOjb8UjdThEsYk15dxSZ030Rw15FppMcFbY7hVH1htfj55K1eS890Sqi5cAEAxN8d+5EicJk/CvFMno713a1OUW0Z+Rin2rpbYONZ9Mev1es7FHePQujUk7N2F/mI9PkVlQ6cBQxj66His7BpvUoncXIUxyfVlXBK8ieaoVRbpbUnU9va0eewxnCZPpmD9BnIWLqTs2DHyVq4kb+VKrAcOxGnyZKwHhsm4uOuwcbS4ZtB2iaIoeAb0wNG9M8knAqgsO4K2/Ah6XREJu1ZzOuZ3/EIHETjiXtx8O9+ClgshhBDG0WLWNm2JFFNT7O+5G++VP9Bh6RJshw0DlYri7dtJefxxzowcSe7KlehkcfZGk5dRCootppZhmNs/jqnVCBR1W3RVVcRu/ZOlr73I92/8g7jtGrRVUqhTCCGuJSUlhUcffRQPDw/MzMzo0KEDzz//PNnZ2YCh5Ndfy3nFx8ejKMpV1SQWLlyIubk5paWlJCUlMW3aNDp27IilpSW+vr7MmjWLioqKW3VqLZoEb7eAoihY9e1L+8/m4rt+HY6THkFlZUV5winS3/wnpyIHkzn3M6qysozWhvTidPam7SW9OP36O7dglyY6ACiKCWrzACwcJjL6lXfxHxiBSm1CWsIJ1n72IV89M5UdPyylKDenaRsthBDNUGJiIv369SMhIYFly5Zx6tQp5s+fz+bNmwkJCSEnJ4fIyEhOnDhBevrle0t0dDSenp5oNJoax4uOjiY4OBhLS0vi4+PR6XR8+eWXHD9+nI8//pj58+fz2muv3eKzbJlkzFsT0RYWkvfjT+R8t5iq1DTAkKmzu+cenKZMxqJr10Z7r1UJq5i9azY6vQ6VomJWyCzu63xfox2/ubnWRIfivFyObF7H4Y1/UHwxaFOp1XQeEEbgiHvx6OJXr65sGZMkjEmuL+NqqWPe0ovTSS5IxsvOCzdrt0ZqYd3uvPNOjh07xsmTJ7G0tLzcjvR0fH19mTRpEh9++CGOjo4sXryYBx54AIAJEybQp08f5syZw5EjR/D29gagQ4cOTJ06lbfeeqvW9/v3v//NF198QWJiorFPrVky2pi3vLw8fv75Z7Zt28bZs2cpKSnBxcWFwMBAoqKiCA0NvamG307Utra0mToFp0cepnDTJnIWLqL00CHyf/6Z/J9/xiokGKfJk7EZNOimCtCmF6dXB24AOr2O2btmE+oReks+/E0hIMwDrwCnWic6WDs4EjL2QfqPup+EvTs5uG4NqSdiObFzKyd2bsW1oy+BI+7FL3QQJlLiRQjRTNzqP8JzcnJYv349c+bMqRG4Abi5uTFx4kRWrFjB559/TlBQENHR0dXBm0ajYcaMGWg0GqKjo5k6dSqJiYkkJycTGRlZ53vm5+fjJFUZ6qVeUUFqaiqPPfYY7u7uvPPOO5SWltK7d2+GDBlC+/btiY6OZtiwYQQEBLBixQpjt7lVUUxMsBsxAu/ly/BevgzbO0eAWk3Jrt2ce+ppEu++h9zly9FdLD3SUMkFydWB2yU6vY6UwpTGaH6zZeNoQbuujnVOdlCbmOAXOogH3/4XD7//Kd0ihqI2NSXjzGnWf/EJXz0zhW3LFlGQlXGLWy6EEDXV9Ue4MYfBJCQkoNfr8ff3r/V5f39/cnNzyczMJDIysrqLNDY2lrKyMgIDAxk0aFD1do1Gg4WFBcHBwbUe79SpU3z22Wc8+eSTxjidVqdembfAwEAmT57M/v37CQgIqHWf0tJSfvnlFz755BNSUlJ46aWXGrWhtwPL3r1p37s3lefPk7NkKXkrV1Jx5gzpb80m8+NPcJgwAceJEzFt61rvY3rZeaFSVDUCOJWiwtPW0xin0CK17ejLiKdfYNDEqRyL3sihDb9TmJXJ3l9WErP6JzoFBRM44h7aB/SQ2cFCiFvuWn+EG7sHpT4jqyIiIpgzZw5paWloNBoGDhyIWq0mPDyc+fPnA4bgLTQ0tNYlK8+fP8+IESO4//77efzxxxv9HFqjemXeYmNj+de//lVn4AZgaWnJgw8+yK5du5g6dWqjNfB2ZNquHW1feZlOGg1tX3sV0/bt0ebnk/3VV5waOpTzL79M6fHj9TqWm7Ubs0JmoVIM/6svpdtba5fpzbCys6f/qHE8Nvd/jPzHa3h264leryNh705+ePs1Fr/8N45sWkelzA4WQtxCl/4Iv5Kx/wjv1KkTiqIQFxdX6/NxcXE4Ojri4uJCWFgYZmZmREdHEx0dTXh4OABBQUFkZWWRmJiIRqNh8ODBVx0nNTWVyMhIQkND+eqrr4x2Pq2NTFhoAfRaLYV//knOokWU7ttfvd0qKAinKZOxiYhAUauveYz04nRSClPwtPWUwK0BspKTOLh+DbHboqkqLwfA3NqagPCh5JhYMGr8AzKgXDQ6mbBgXC1xwkJTTDyLiori+PHjJCQk1Dlh4YsvvgBg0KBB+Pn58csvv7BmzRr69+8PwJAhQwgODubdd99lx44dNcbGnz9/nsjISPr27cuSJUtQX+c+1toZfYWF1NRUtm/fTkZGBjpdzVTuc88919DD3TItNXi7UunRY+QsWkTBunVQVQWAqZcXTpMm4TBmNCrrxl8KSkBZURHHt2zi4Po15F+4PM7EO7Affe8cSYeegdKlKhqNBG/G1RKDN7j1f4QnJCQQGhqKv78/77zzDh07duT48ePMmDGD8vJydu/eXT3BYNasWXz88ceAYbKDiYlhVNbbb7/Nhx9+iE6nIzc3t/p6Pn/+PBEREXTo0IFFixbVCNzc3G7PBINRg7eFCxfy5JNPYmZmRps2bWrcsBRFadZTfFtD8HZJZXo6uUuXkrviB3QFBQCo7OxwuH8cTg8/jKm7exO3sHXS63ScObSf/Wt/Jfnowertjh7t6T38brqFD8HcyqoJWyhaAwnejKulBm9N4ezZs8yaNYt169aRk5ODm5sbo0ePZtasWbRp06Z6P41GQ2RkJCNGjOCPP/6o3r5lyxYiIiKIiopi3bp11dsXLlxY5xCr26hDsAajBm+enp489dRTvPrqq6huooRFU2hNwdsluuJi8n75hZzFi6k8m2zYqFZjFxWF05TJWPbs2bQNbKUqKyv5Zfn3OFaUELftTyouzgY2tbCkW/hgeg+/hzbtZVKIuDESvBmXBG+iOTLq2qYlJSU88MADLS5wa61U1tY4TZyI44MPUqTZQs6iRZTs2UPB2rUUrF2LZZ8+OE2ejO3QIdcdFycaxszOgfC7HqLPXQ9weONGEvdvJi/9PIfW/86h9b/j1aM3gVH34NM3CJVKfvdCCCEaR4ODt2nTprFy5UpmzpxpjPaIG6SoVNgOjsR2cCRlcXHkLFxE/tq1lB44wPkDBzBt1w6nSY9gP3YcahsZF9dY4nels21ZAnq9MzCeoDGQk7KLxP17ST56iOSjh7BzcaXXsLvoMXg4lratI+MrhBCi6TS421Sr1XLPPfdQWlpKjx49rkrp/+c//2nUBjam1thtei2VGRnkfv89ectXoM3LA0Bla4vjhPE4PvIIpm3bNm0DW7DKykp+/ekPLmyx4cpPkKKCSXNC0Vbmc3jTHxzdvJ6yokIATEzN6Bo2iMAR99K2o28TtVy0BNJtalzSbSqaI6N2m7733nusX7+erhfX3vzrhAXRfJi6uuL6wgs4P/kk+at/JWfRIirOnCH7fwvIXrgI+7vvxunRqY26jurtpKpExV//9NHrID+jlHZd2zLooSmEjHuQEzu2cnDdGjKSTnNcs4njmk14dPGn94h76DIgFLWJ3JyFEELUX4ODt48++ohvvvmGKVOmGKE5whhUlpY4PjABh/H3G8bFffMNJfv2kb96NfmrV2MdGorTo49iHRYqAXgDmFjpUBSuyrzZu16uh2RqZk73yGF0ixhK6sl4Dq1fw8nd20k9GUfqyTi2ODjSc+gIeg4ZgY1Tm1reRQghhKipwbMOzM3NCQsLM0ZbhJFdGhfXYcl3eK/8Abu77gSViuKdO0l57DHOjBpN3i+/oK+oaOqmtggmlnrueLAzlwqfKyqImOhX63qqiqLQrqs/dz83g8f/+y0h4x7C2sGR4rxcdv24jK+ffZQ1n/6L8/Gxt+00eSGEEPXT4DFv7733HmlpacydO9dYbTKa223MW31UnDtHzuLF5P34E/qSEgBMXF1xfORhHCdMQC2/p1pdOSapvEhLfkYp9q6WtQZuddFWVZKwZycH1/9O6onY6u0u3j70Hn43/gPDMTWXMTK3IxnzZlwy5k00R0at8zZmzBj+/PNP2rRpQ7du3a76Ylm1alXDW3yLSPBWN21+Prk//EDu4u+oyswEQGVlhf24sThNmoxZ+3ZN3MLmpbFvrhfOnObQ+jXEb99CVaUh82lubU238KH0Hn4Xju7y+7+dSPBmXBK8iebIqBMWHBwcuO8+466nJm49tb09zo8/TpvJk8n/fS05335L+cmT5C7+jtwlS7EbEYXT1Eex7NH9usdKL04nuSAZLzsvWUe1ntp29CXqqecZNHEqxzSbOLxxLfkX0jmwdjUH1q6mQ89Aekfdg0+fflIzTgghbnMNDt6+/fZbY7RDNBOKmRkOY0ZjP3oUxTt2kvPNNxTv3EnB2j8oWPsHVkFBOE2dik1EOEothZqbYvHk1kRbZYZ7l0j8QkeQmXycQ+t/58yh/Zw9cpCzRw5i5+JKz6F30mPwcKzs7Ju6uUKI28ilJbByc3NxcHC44eOUlJTwyCOPsHHjRgoLC8nNzaV379688MILvPDCC43W3tZMlkkQtVIUBZuBYXh9s4COv/yM/ahRYGJCSUwM5555hsS77yH3hx/QlZdXvya9OL06cAPQ6XXM3jWb9OL0ut5GXCF2RyqLX9vJ6o8P8t0buykracd9M99i2qdf0+/e+7CwsaUgM4Ptyxbx1dOTWTvvI1JPxssEByFEo4uIiLgqkAoNDSUtLQ17+5v7w3HRokVs27aNnTt3NsrxblRKSgqPPvooHh4emJmZ0aFDB55//nmys7Or95k5cyZ+fn41XhcfH4+iKFdV3Vi4cCHm5uaUlpaSlJTEtGnT6NixI5aWlvj6+jJr1iwqGmlCYL2CtxEjRrB79+7r7ldYWMgHH3zAf//735tumGg+LPz88PjgfTpt2kibx6ahsrWl4swZ0v85i1ODh5D5+edU5eaSXJBcHbhdotPrSClMaaKWtxxFuWVolsRXlx3R60GzNJ6i3DIc2roR/vCjPPHFQqKefoG2Pp3RVlURty2aZW++xJJXX+Bo9AYqK8qv/SZCCHETzMzMcHNzu+mSUqdPn8bf35/u3bs3yvFuRGJiIv369SMhIYFly5Zx6tQp5s+fz+bNmwkJCSEnJweAyMhITpw4QXr65SREdHQ0np6eaDSaGseMjo4mODgYS0tL4uPj0el0fPnllxw/fpyPP/6Y+fPn89prrzVK++sVvN1///2MHTuWgIAAXnnlFVauXMmOHTvYv38/mzZtYu7cuYwfPx53d3cOHDjAvffe2yiNE82LqZsbri+9RKfoaNq+OhMTD3e02dlkzf2MU5GDcf7vT3jk1vwQqhQVnrY3vkB7enE6e9P2tvrsXV5GaZ0Ffy8xNTOne8RQHn7vYx6a8xEBgwajNjUl48xpNsyfy1dPTUaz+H/kpp2/xa0XQrQmU6ZMYcuWLXz66acoioKiKCQlJaHRaFAUhbyLK/YsXLgQBwcH1qxZQ9euXbGysmLcuHGUlJSwaNEivL29cXR05LnnnkOr1QKGjN5HH33E1q1bURSFiIiIWtuQnJzMqFGjsLGxwc7OjvHjx3PhwgUA8vPzUavV7Nu3DwCdToeTkxPBwcHVr1+yZAmennXfe6ZPn46ZmRkbNmwgPDwcLy8v7rzzTjZt2sT58+d5/fXXARg4cCCmpqY1AjWNRsP06dPJyckhKSmpxvbIyEjAkPT69ttvGT58OD4+PowcOZKXXnqp0SZ11mvM27Rp03j44YdZuXIlK1as4KuvviI/Px8wdK8FBAQQFRVFTEwM/v7+jdIw0XypbaxxmjwZx4kTKVi/npxvvqXs+HHKf/yVjxWFmM7w2wAVCZ5qZoXMuuFJC7fT+DkHV8vrFvy9knunrrh36kr4I9M4Fr2Rwxv/oCDzAvt//4X9v/+CV4/e9B52F779BqBSywQHIZoLvV6PvrT0+jsagWJpWa8s16effsrJkyfp3r07b7/9NgAuLi41ApVLSkpKmDt3LsuXL6ewsJD77ruPMWPG4ODgwNq1a0lMTGTs2LGEhYUxYcIEVq1axcyZMzl27BirVq3CzMzsqmPqdLrqwG3Lli1UVVUxffp0JkyYgEajwd7ent69e6PRaOjXrx9Hjx5FURQOHjxIUVFR9evCw8NrPb+cnBzWr1/PnDlzsLSs+R3r5ubGxIkTWbFiBZ9//jnW1tYEBQURHR3NAw88ABiCtBkzZqDRaIiOjmbq1KkkJiaSnJxcHbzVJj8/Hycnp+v+/uuj3hMWzM3Nefjhh3n44YerG1FaWkqbNm1kKvttSjExwf7uu7G76y5KYmLI+eZbijQa+p+E/ie1qHt0xa2tDXofLUoDA4i6xs+FeoS2yhmsNo4WRDzsh2ZpPHrdtQv+XsnKzp7+o8bR794xJB06wOGNa0k8uI/ko4dIPnoIG0cnegyJoseQKGydnG/R2Qgh6qIvLeVEn75N8t5dD+xHsbK67n729vaYmZlhZWWFm9u1v28rKyv54osv8PU1rNc8btw4vvvuOy5cuICNjQ0BAQFERkYSHR3NhAkTcHJywsrKqroLtjabN2/m6NGjnDlzpjp7tnjxYrp160ZMTAxBQUFERESg0Wh46aWX0Gg0DBs2jPj4eLZv386IESPQaDS8/PLLtR4/ISEBvV5fZ7LJ39+f3NxcMjMzcXV1JTIykpUrVwIQGxtLWVkZgYGBDBo0CI1Gw9SpU9FoNFhYWNTI/l3p1KlTfPbZZ3z44YfX/H3WV4Nnm15ib2/fZIMMRfOiKArW/ftj3b8/5adPk7NwIfmrf0V7NJbzzz+PqZcXTpMn4TBmDKp6fHEA1xw/1xqDN4CAMA+8ApxuqOCvSqXGtWMPgkZ3YsB9BZyOieZo9EaKcnPY9eMydq9agW/fAfQafhcduveqdaawEEI0lJWVVXXgBtC2bVu8vb2xsbGpsS0jI6Pex4yLi8PT07NGt2dAQAAODg7ExcURFBREeHg4CxYsQKvVsmXLFoYPH46bmxsajYaePXty6tSpOrtkL6nvZK+IiAjmzJlDWloaGo2GgQMHolarCQ8PZ/78+YAhGxcaGoq5uflVrz9//jwjRozg/vvv5/HHH6/37+Fabjh4E6I25r6+uP/f/+Hy/PPkfv89uUu/pzI5mQv/9w5Zcz/D4aEHcZo4ERPna2eBvOy8UCmqGgHczY6fawlsHC0aFLRdErsjtXrCg6JAxMPDeeLziZzau5PDG//gXNwxTsXs4lTMLhzc3Ok19E66RQzF0laKVQtxKymWlnQ9sL/J3rux/bXnTVGUWrfpdDX/GL9ZgwYNorCwkAMHDrB161beffdd3NzceP/99+nVqxceHh507ty51td26tQJRVGIi4tjzJgxVz0fFxeHo6MjLi4uAISFhWFmZkZ0dDTR0dHV3bFBQUFkZWWRmJiIRqPhySefvOpYqampREZGEhoayldffdVo5y9/fgujMHF2xuW55+gU/Sdt//kmpl5eaPPzyf5iPqcGDyHtzTcpP326zte7WbsxK2QWqosLh14a89Zas243o66ZqmVFWvzCwpnw1vtM/vC/9I66BzNLK/LS09iy5Bu+vFhu5FzcMSk3IsQtoigKKiurJnk0ZFanmZlZ9SSDW83f35+UlBRSUi5XKoiNjSUvL4+AgADAsGBAz549mTdvHqampvj5+TFo0CAOHjzImjVr6hzvBtCmTRuGDRvG559/Tulfxh+mp6ezdOlSJkyYUP37srS0ZMCAAWg0GrZs2VKd0TM1NSU4OJgFCxaQkpJy1Xi38+fPExERQd++ffn2229RNWKPhwRvDZSWX8rO01mk5TfNgNOWRmVlhdNDD+H7x1razf0Uy1690FdUkLfyRxLvvoeUp56meO/eWoOH+zrfx/qx6/km6hvWj13faicr3Kz6zFR19uzAkEef4sn5ixj2xLO4evuirawkbls0K96aycK/P83+33+htLDgFrdeCNEceXt7s2fPHpKSksjKymr0zNm1DB06lB49ejBx4kQOHDjA3r17mTRpEuHh4fTr1696v4iICJYuXVodqDk5OeHv78+KFSuuGbwBzJs3j/LycqKioti6dSspKSmsW7eOYcOG0a5dO+bMmVNj/8jISJYvX05ZWRl9+vSp3h4eHs5nn31WPbHhkkuBm5eXFx9++CGZmZmkp6fXKDlyMyR4a4AVMcmEvf8nD329h7D3/2RFTHJTN6nFUNRq7IYPx3vFcjp8vxSboUNAUSjSaEieNJmk+8dTsHYt+qqqGq9zs3YjyC1IMm7XcGmm6pXqmqlqZmFJzyEjePj9T3hozkf0GDwcU3MLclLPoVn8P758ahK/z/03KbFHJRsnxG3spZdeQq1WExAQgIuLC8nJt+5+pygKq1evxtHRkUGDBjF06FB8fHxYsWJFjf3Cw8PRarU1xrZFRERcta02nTt3Zt++ffj4+DB+/Hh8fX154okniIyMZNeuXVfNCo2MjKSwsJCwsDBMTC6POAsPD6ewsLC6pMglGzdu5NSpU2zevJn27dvj7u5e/WgMDV6YHiAvL48ff/yR06dPM2PGDJycnDhw4ABt27alXbvmu4D2zSxGnJZfStj7f6K74relVhS2z4zE3b7xxxHcDsrPnCFn0SLyf/4F/cWVGkw9PHCaMhmHsWNRWVs3cQvr1twWDo/dkXrVTNWAMI96vba8pIT4HVs4smkdGUmXu7Id3dvRY0gU3cKHyFJct1hzu75aG1mYXjRHRl2Y/siRIwwdOhR7e3uSkpJ4/PHHcXJyYtWqVSQnJ7N48eIbbnhzdiaruEbgBqDV60nKKpHg7QaZd+yI+1tv4fLcc+R+v4zcpUupTE3lwrvvkTnvvzg+8ACOD0/E1NW1qZva7N3MTFVzKyt6DbuTXsPu5ELiKY5sWkfcji3kpp1n65Jv2L5sMZ37h9Bz6Ag8A3rITFUhhGhiDf4W/vvf/86UKVNISEioERneddddbN26tVEb15x0dLZGpYAtJThiGBekVhS8netX+kLUzcTJCZdnp9Mp+k/c3noLsw4d0BUUkP3VV5waMpTUV1+j7OTJpm5ms2fjaEG7ro71CtyKcss4dyKXotyyGtvb+nRi2BPP8tTFsXFtfTqj01ZxYtc2Vv7f63zz4pPsXf0jJfl5RjoLIYQQ19PgzFtMTAxffvnlVdvbtWvXaAPxmiN3e0veu68Haavf4kn1b/ygjcRp2N8l69aIVBYWOD4wAYfx91MUHU32N99Sun8/+T//TP7PP2M9cCBtHp2KVUhIk6yF11pcXVbk6i5WM0sreg4ZQc8hI7iQeIqjf64nbruGvPQ0tn2/kB0rltApKJieQ0bg1b2nZOOEEOIWanDwZm5uTkHB1TPSTp48WV0TpbWaEORF2bEcLM5WMNlkPWg2QfY4CHsB2gY0dfNaDUWlwnbIEGyHDKH00CGyv11I4caNFG/fTvH27Zj7+dHm0anY3XknSh3jgdKL00kuSMbLzksmO1yhrrIiXgFOdWbs2vp0oq1PJwY9/Cgndm7jyOZ1pJ86ycnd2zm5ezv2bd3oMdgwNs7GsXGWfhFCCFG3Bv+5PHLkSN5++20qKysBw6yQ5ORkXnnlFcaOHdvoDWxuLKb8DJNWg08E6LVwZAV8EQLfT4Czu5q6ea2OZe/etP/0E3zXr8Px4YdRLC0pj48n9eVXODV0GNkLFqAtLKzxmlUJq4j6KYppG6YR9VMUqxIaZyHg1qA+ZUXqYmZhSY/Bw5k45z888sFcekfdjZmlFfkX0tm+bBFfPTOFnz+YTcKenWirKo10BkIIIRo82zQ/P59x48axb98+CgsL8fDwID09nZCQENauXYt1M54h2OgzjFIPwvZPIHY1cPHX6BkMA1+EzsNBupIanTYvj9zlK8hZugRtZhZgqCXncP/9OE16hGwHNVE/RV21MsP6sesbNQPXUmcDFuWWsfi1nTUCOEUFk+aE3tDKDpVlZZzYvZ0jm9eRdjK+erulrR0BgyLpFjEMFy/vRmj57aWlXl8thcw2Fc1RQ66VGyoVArB9+3aOHDlCUVERffr0YejQoTfU2FvJaB/YrFOwcy4cXgbaCsM2F38Y+AJ0Hwtq+fJtbLqKCgp+W0POwm8pTzhl2KhWUxnejze99pHoXnNM3DdR3xDkFlTLkW5MS7653kxZkWvJPp/C8S2bid2ymeK83OrtbX060z1iKH5h4Vhcsd6hqFtLvr5aAgneRHN0S4K3lsjoH9jCdNj9OcR8AxUXu/LsPSHkWejzCJg136xkS6XX6ynevp3sb76hZNfu6u3HvWBNfxUHOikoKrVk3v6iKLfshsqK1IdOqyXp8AGORW/k9P496C4usaM2NaVz/1C6RQylQ/deMsnhGlr69dXcSfAmmiOj1nmbO3durdsVRcHCwoJOnToxaNAg1Gp1Qw/d8tm6wbC3YeDfYd8C2P0F5KfAuldgywcw4Cno/zhYyaDuxqIoCjZ33IHNHXdQFhdH9rffkv/773RL1tEtWcd5J1A9cA+uKikyeyUbR4tGD9ouUanV+PQJwqdPECUF+cRt03BMs5Gs5CTid2whfscWbJ1d6BY+hG7hQ3FoKxNKhBCiIRqceevYsSOZmZmUlJTg6OgIQG5uLlZWVtjY2JCRkYGPjw/R0dF4enoapdE36lb+tQVAZSkc+t7QpZqbZNhmagV9p0DIdLBvb/w23IYq09NJ+WY+ZT/9hlJcAoDawQGHBx/A6aGHMGmEWdGSGampKLeMvIxSHOrI5On1ejLOnOZo9Ebid2goLy6ufs4zoAcBgwbTeUAY5lZSNxHk+jI2ybyJ5qgh10qD+y3effddgoKCSEhIIDs7m+zsbE6ePMmAAQP49NNPSU5Oxs3NjRdffPGGT6DVMLWEoGnw7H4Y9w249YDKEkPX6qe94OenISP++scRDWLq5obPa2/RdctW2r72Kqbt2qHNyyP7i/mcGjyE1Ndel6K/jSh2RyqLX9vJ6o8Psvi1ncTuSL1qH0VRaOvTiaHTnuap+d9x93Mz6NAzEBSFlNijrJ//KfOfeJg1n/6LxIMx1V2tQoimlZKSwqOPPoqHhwdmZmZ06NCB559/nuzsbABmzpyJn59fjdfEx8ejKApTpkypsX3hwoWYm5tTWmqY3T5y5Ei8vLywsLDA3d2dRx55hNTUq78/xNUanHnz9fXlp59+onfv3jW2Hzx4kLFjx5KYmMjOnTsZO3YsaWlpjdnWm3bLM29/pdfD6c2GGapJ2y5v73q3YXKDZ/9b36bbgL6qisJNm8n59ltKDx+u3m4dFobT1KlYh4U2uOivZEYMbnb2akFWBnHbNMRu/ZOc1HPV263sHfALCyfgjkhcO/redkWZ5foyLsm81U9iYiIhISF06dKFd955h44dO3L8+HFmzJhBRUUFu3fvJiYmhhEjRpCWloabm2EIxBdffMF7772HSqUiKSmp+niTJ08mKSmJLVu2APDxxx8TEhKCu7s758+f56WXXgJg586dt/xcmwOjZt7S0tKoqqq6antVVVX1CgseHh4U/qX21s1ISkpi2rRpdOzYEUtLS3x9fZk1axYVFRWN9h63hKJAp6EwZQ08thn87gEUOPE7LBgG394FJzdwVSEucVMUExPsRkThvWI5HZZ9j+1wQxmX4h07SHnsMc6MHEXeT6vQtbTrqRm4mbpxAHbOrgwYM54p//mCiXP+Q+CIe7G0taMkP48Da1ez5NUXWPTSdPau/pHC7CwjnIEQLUddy9oZy/Tp0zEzM2PDhg2Eh4fj5eXFnXfeyaZNmzh//jyvv/46AwcOxNTUFI1GU/06jUbD9OnTycnJqRG8aTQaIiMjq39+8cUXCQ4OpkOHDoSGhjJz5kx2795dXUdW1K3BwVtkZCRPPvkkBw8erN528OBBnn76aQYPHgzA0aNH6dixY6M1Mj4+Hp1Ox5dffsnx48f5+OOPmT9/Pq+99lqjvcct174fPLAUpu+FwIdBZQpnd8D398P8gXBkJWivDpLFzbEKDKT93E8NRX8nPYJiZUV5QgJpr7/OqcFDyPriC6pyc69/IAGAg6slf02KKSqwd23YsnGKouDWqQuDpz7Jk/MXM/rlf9IleCBqU1OyzyWz7fuFfDV9Kiv/7zWOaTZRXlLSiGchRPNXn+EJjSknJ4f169fzzDPPYGlZ8/Ps5ubGxIkTWbFiBVZWVgQFBREdHV39vEajYciQIYSFhVVvT0xMJDk5uUbw9tf3W7p0KaGhoZJtrocGB28LFizAycmJvn37Ym5ujrm5Of369cPJyYkFCxYAYGNjw0cffdRojRwxYgTffvstw4cPx8fHh5EjR/LSSy+xalUrqJzv0gVG/ReeP2woKWJmAxeOwarH4LNA2Ps1VMiNqrGZeXri9tprdNZE4zrjJUzatkWblUXmp3M5FTmYtLfeojzxTFM3s9mzcbQg4mE/lIvfJJfqxtV3JmttmQS1iQm+fftz74szeerL7xj2xN9o798d9HqSjx1h/Ref8MUTE/n1o3c5uXs7lRXlxjg1IZqNupa1M2YGLiEhAb1ej7+/f63P+/v7k5ubS2ZmJpGRkdWZt9jYWMrKyggMDGTQoEHV2zUaDRYWFgQHB9c4ziuvvIK1tTVt2rQhOTmZ1atXG+2cWpMGlwpxc3Nj48aNxMfHc/LioO+uXbvStWvX6n3qiqwbU35+Pk5O1y65UV5eTnn55S/2S2uyVlZWNr+0rJUrDH4LQp5Htf8bVDFfoeQlw9qX0GveRxf0OLq+08DSoYkb2spYWmI3aRK2Dz5I0YYN5C1aTHlcHHnLV5C3fAVWEeE4TJqEZb9+NcZdXbp+KisruVBygeTCZLxsvWhr1bapzqTJdO7vgntnOwqySrFztsTG0bxen6/4XelsW5aAXm8YUXDHg53xC6lZNkRtZo7/oMH4DxpMfsYFTuzcyomdW8hNPU/C3p0k7N2JqYUlvv0G0CVkIJ7deqE2afDXWrNz5fUlGl9L+71ea3iCsUr+VL9PPYbxREREMGfOHNLS0tBoNAwcOBC1Wk14eDjz588HDMFbaGgo5ubmNV47Y8YMpk2bxtmzZ5k9ezaTJk1izZo1t90414a64W85Pz+/q2aY3CqnTp3is88+48MPP7zmfu+99x6zZ8++avuGDRuwatYlCfxRd34fr+xt+GasxbokC/WW99Bv+5gk50hOu0RRZia14hqdosDkSVieOYPj1m3YxMVRotlCiWYLZe3akXvHQAp79oQrahi+t+Y9VpeuRo8eBYVRlqPoZ96vCU+iiSXUb7eqUoV0jTVg+ILW62Hr9yc5mXoQE8tr3CzMrHEKvxObvBwKk05RdDaRypIi4rdriN+uQWVugY1nR2y9fbFwcWvxN4CNGzc2dRNapZIW1u1+aXjCXycGNXR4QkN06tQJRVGIi4tjzJgxVz0fFxeHo6MjLi4uhIWFYWZmRnR0NNHR0YSHhwMQFBREVlYWiYmJaDQannzyyauO4+zsjLOzM126dMHf3x9PT092795NSEiI0c6tNWjwbFOtVsvChQvZvHkzGRkZ6HS6Gs//+eef9T7WzJkz+eCDD665T1xcXI0g8fz584SHhxMREcH//ve/a762tsybp6cnWVlZTTPb9EboqlBif0G9ay5KRiwAepUp+h7j0QY/C86dm7iBrVfFmTPkLVlC4epf0V+8jtSurjhMfAir0aNZtXM9HxV8hI6a66j+Pur32zID1xCpJ/NY89nRq7bf81wPPDo71Pm6otxyCjJLsXMxZPj0Oh1pp05wctc2EvbsoPRidh3AxqkNnfqH0rl/KG6durSoFR0qKyvZuHEjw4YNk/E/RlBQUICzs3OLmm1qrGXtriUqKorjx4+TkJBQY9xbeno6vr6+TJo0iS+++AKAQYMG4efnxy+//MKaNWvo399QPWHIkCEEBwfz7rvvsmPHDkJDQ+t8v+TkZDp06EB0dDQRERFGPbfmyKjLYz377LMsXLiQu+++G3d396v+sv3444/rfazMzMzqWjF18fHxwczMDIDU1FQiIiIIDg5m4cKFqBr4ZdzkpUJuhl4PCRth+8eQfGkatQJ+dxtWdGjft0mb15pV5eaSt3w5OUu/R5tlmPGoWFpyprcv/+4dR6aDcddRbY1upMRI7I7U6nE/igIRD9e8eem0WpKPHSZ+x1YS9u6kovRydsXGqQ1dBoTRJXggHl38mn0gJ6VCjKullgox5rJ2tUlISCA0NBR/f/+rSoWUl5eze/fu6uFLs2bNqr7/5+TkYHJx+MLbb7/Nhx9+iE6nIzc3t/p63rNnDzExMQwcOBBHR0dOnz7Nm2++yYULFzh+/PhV3au3A6MGb87OzixevJi77rrrphrZUOfPnycyMpK+ffuyZMmSG1p+q0UHb1dK3gM7PoETay9v877DUCvOdwhXTf8D0vJLOZNVTEdna9ztjZdqb810FRUU/LaGnIULKU8w9A/qFNjbRWFNfxUn24HKCOuotlYNySQ0NNjLyyjgxM69pJ/eT/LRfVSUXi5dYuPUhs4DQukSPJB2XfybZSAnwZtxtdTgrSmcPXuWWbNmsW7dOnJycnBzc2P06NHMmjWLNm3aVO93qQzIiBEj+OOPP6q3b9myhYiICKKioli3bl319qNHj/L8889z+PBhiouLcXd3Z8SIEbzxxhu0a9fulp5jc2HU4M3DwwONRkOXLl1uqpENcf78eSIiIujQoQOLFi2qEbhdKgpYH60meLskIx52fApHfwDdxbIibj1g4IvgPwrUhr98VsQk8+qqo+j0oFLgvft6MCHIqwkb3rLp9XoKtm4l4cOPsE64PMjrlLuC5UP3M2TKGyhyw62X+mYSzp3IZfXHB6/aPvrFQNp1dayx7a8Zujse8MXSOp2Tu7dzet+eGhk5a0cnOvULplP/EDwDejSbyQ4SvBmXBG+iOTJq8PbRRx+RmJjIvHnzbtlg4IULFzJ16tRan2tI81td8HZJXgrs+i8cWGRYfgvA0RtCnyPN5z7CPtyJ7opfk1pR2D4zUjJwN+HSzXVI585cWPQ/ytduQrk4g83E1RXHiRNxGH8/Jo6O1zmSqI/6Zt6ut19VZSVnjxzg5O4dnN63h/KSy2usmltb49OnP52DQvDu1QfTJrzRSvBmXBK8iebIqMHbmDFjiI6OxsnJiW7dul31xdKca6+12uDtkpIc2PsV7PkSSnMAqLBw5pPCISzRDqUA6+pdlz0eTIhvm7qOJK7jrzfXquxsclesIPf7ZZfHxVlYYD9yJE6THsG8U6cmbnHLV59u1oZk6KoqK0k5dpiEmF2c3reHkvy86udMTM3o0CuQTv2C8ekThJW9gzFOqU4SvBmXBG+iOWrItdLgPgIHB4dapw2LZsDKCSJmQujf4MB3sPMzzArO8bLpCp42+ZWl2iF8U3Un2YoT3s7NuVRKy2PSpg0uzzxDm8ceo2DtWnIu1Yv74QfyfvgB64EDcZo8CeuwsGuOsUovTie5IBkvOy8ZN/cXAWEeeAU4XbObtSElFcqKtJha+RJ6f3eGPvYMqSfjObV3F6didpGfcYHT+/Zwet8eUBTcO3fFt09/fPsNoE17rxZfgkQI0bI1OPPWkrX6zNtfaSvh6I/kb/oQ+6JTAJTrTTjnNQrfUa+Bs2SDbtT1MiN6vZ7SffvIXrSIos1/VkcTZj4+OE16BPtRo1D9ZcmZVQmrmL1rNjq9DpWiYlbILO7rfN8tOZ/WpD4ZumvNXNXr9WSePcOpmF2c3reXjKTTNV5r79oWn7798e0zgHb+3TAxQmZMMm/GJZk30RwZtdu0JbvtgrdLdDpyDv2Gya5Pscvcf3GjAv73GmaotpMyIw3VkJtrRUoKuUuWkPfjT+iKDWOsVPb2OI4fj+PEhzB1cyO9OJ2on6LQ6WvWjJOZqzfmWhMhGjpztTA7i8QDezm9bw/Jx4+gvaI6v6m5BV49etGxdz86BvbFztm1UdovwZtxSfAmmiOjdpsC/Pjjj/zwww8kJydTUVFR47kDBw7cyCGFMalUOPUZBX1GwdldhjIjJ9dB3K+GR8dBEPYC+A6utcyIuDlmnp60ffVVnP/2N/J/+omc75ZQee4c2V9/TfY332AXFUXa3f1qBG4AOr2OlMIUCd5ugI2jRZ2zV+u71FBRbhl5GaU4uNrQa9hd9Bp2FxVlpZw9eojE/Xs5c3AfxXm5l7tXgTbtvegY2A/vnn1o5xeAycUalUII0ZgaHLzNnTuX119/nSlTprB69WqmTp3K6dOniYmJYfr06cZoo2hMHUIMjwuxF8uMrIQzWw0Pt56GTNwVZUZE41Hb2OA0eTKODz9MUXQ0OYsWUxITQ8HatVivXcs77RR+D1LY01VBp1JQKSo8bT2butmtTn3GxdXVrWpmYUnnoBA6B4Wg1+nIOHuGpEP7STy4j7ST8WSfSyb7XDL7fluFiakZ7QO606FnIN49A2nj2UHGygkhGkWDu039/PyYNWsWDz74ILa2thw+fBgfHx/++c9/kpOTw7x584zV1pt223abXkte8sUyI4uvKjNC74lgKmn+2jRWt1bp8ePkLv6O/LVr4WJ3XJYtbOirpt9jrzCq3yON1WRxhWuNi6tPt+rlrNzlbtnSokLOHjlI0qEDnD1ygKLcnBrvaWXvgGdADzy79aB9QA+cPNrXGcxJt6lxSbepaI6MOubNysqKuLg4OnTogKurKxs3bqRXr14kJCQQHBx83eWumpIEb9dQnG0oM7L3SyjNNWyzdoXgp6DfNLB0aNLmNTeNfXOtyswkd9lyspd9jz43DwDF3By7e+7G6eGHsfD3v+n3EDXVNS7ueuVGrjXZ4VJQZ+9iQXlxhiGYO3KQc7HHqKoor3G8awVzErwZlwRvojky6pg3Nzc3cnJy6NChA15eXuzevZtevXpx5syZBhXMFc2MdRuIfBXCnqsuM0LBOdj8Nmz7GPpNheBnwM69qVvaKpm4uODy3N9o8+QTFKz9g5zvFlMeG0f+T6vI/2kVVv364fjII9gOGYxSxyoAUmakYeoaF3etbtWi3LLqwA0M+2iWxuMV4ERybM5VQV3fu0fT9+7R5GUUkHjwOIWZp7iQGEfqyThK8vM4sWsbJ3ZtAwzBXPuAHngG9MC9i598n4pm59ISWLm5uTg4ONzwcUpKSnjkkUfYuHEjhYWF5Obm0rt3b1544QVeeOGFRmtva9bgRf0GDx7Mr7/+CsDUqVN58cUXGTZsGBMmTJD6b62BmbUh2/b8IRjzJbj4Q0Uh7JwLn/aE1c9CVsJ1DyNujMrcHIcxo+n40090+H4ptneOALWakn37OP/885waNpysr7+mKje3xutWJawi6qcopm2YRtRPUaxKaL7Fsps7G0cLIh72Q7n47XipW9XG0aLOyQ7pp/NrDeqKcsuI3ZHK97P2sevnUo7vaEf3Ic/y7DcrmDDrfULvn4hnt56YmJpRkp/HyV3b2Lzgc5a88hxJPy/lj88+ZM8vqzm25TAF2cUIcatERERcFUiFhoaSlpaGvb39TR170aJFbNu2jZ07dzbK8W5USkoKjz76KB4eHpiZmdGhQweef/75Gj2IM2fOxM/Pr8br4uPjURSFKVOm1Ni+cOFCzM3NKb24lvLIkSPx8vLCwsICd3d3HnnkEVJTUxul7Q3OvH311VfodIZZcdOnT6dNmzbs3LmTkSNH8uSTTzZKo0QzoDaFXg9Aj/GQsMEwQzV5Fxz8Dg4uAf97IOxFaH91mZG0/FLOZBXT0dlaluC6QYqiYNWnD1Z9+lCZnk7u8uXkrfiBqrQ0Mj/6D1nz/ovdvffg9Mgj5LW3r64PB4ZZqrN3zSbUI1QycDeoroLAdWXl9NCgoM4rIJT2Ad1pH9CdEB6kqrKS9FMnSIk9yrnYo5w/EY+2rJSEPTtI2LPj4hFNcWrXEd9+PXHv1AX3Tl2xcZJVUsStY2Zm1qD1xOty+vRp/P396d69eyO06sYkJiYSEhJCly5dWLZsGR07duT48ePMmDGDP/74g927d+Pk5ERkZCQffPAB6enp1eceHR2Np6cnGo2mxjGjo6MJDg7G8mINz8jISF577TXc3d05f/48L730EuPGjWPnzp033f4GZ95UKhUmV3TbPPDAA8ydO5e//e1vmMm0+NZHpYKuI+DRdfDoeuhyJ6CHuN/gf4Nh4T1walP1nWtFTDJh7//JQ1/vIez9P1kRk9y07W8FTN3ccH3hBTpponF/913MA/zRl5eT/+NPnBk1msxHn6ZffBWqKxawvVRmRNw4G0cL2nV1rNG1WldWzt3X/qoqO9cK6vIzSmtsMzE1pb1/d0LGPsj9b77Lk18twe2OezGxDEVl4gmYApXknD9JzOof+fWjd/ny6cl89cxUfvvPe8T8topzcceoLC9r9N+DaFx6vZ7KsrImedS3K37KlCls2bKFTz/9FEVRUBSFpKQkNBoNiqKQl5cHGDJNDg4OrFmzhq5du2JlZcW4ceMoKSlh0aJFeHt74+joyHPPPYdWqwUMGb2PPvqIrVu3oigKERERtbYhOTmZUaNGYWNjg52dHePHj+fChQsA5Ofno1ar2bdvHwA6nQ4nJyeCg4OrX79kyRI8PeuerT99+nTMzMzYsGED4eHheHl5ceedd7Jp0ybOnz/P66+/DsDAgQMxNTWtEahpNBqmT59OTk4OSUlJNbZHRkZW//ziiy8SHBxMhw4dCA0NZebMmezevZvKK2pF3qgbqgeRl5fH3r17ycjIqM7CXTJp0qSbbpRopryC4aHlkBF3ucxI0jbDw60HuX2m88bP1uj0agB0enht1TEGdXGRDFwjUJmb43DfGOzHjKb04EFyvvuOwg0bMTkcz0uHIdMONvRR8WcvhWJrtZQZMZK6snIRD/tdNYP1UlBXn+W6rmRiaoqJVTtMLDqDRTB6vQ69LgddVRpe/pUUZCaRnZJMYXYmhdmZnLyYnVNUKly8OtLWxxcXbx9cO/jg0sEbM0tZDq+5qCovZ+7kcU3y3s8t+hHTekya+PTTTzl58iTdu3fn7bffBsDFxaVGoHJJSUkJc+fOZfny5RQWFnLfffcxZswYHBwcWLt2LYmJiYwdO5awsDAmTJjAqlWrmDlzJseOHWPVqlW1Jn10Ol114LZlyxaqqqqYPn06EyZMQKPRYG9vT+/evdFoNPTr14+jR4+iKAoHDx6kqKio+nXh4eG1nl9OTg7r169nzpw51VmyS9zc3Jg4cSIrVqzg888/x9ramqCgIKKjo3nggQcAQ5A2Y8YMNBoN0dHRTJ06lcTERJKTk2sEb399z6VLlxIaGtook5AaHLz99ttvTJw4kaKiIuzs7GpMdVcURYK324GrP4yZD5Gvw+7PYf9CSD+K49qn2GTqytfau1mpDaccM7R6PUlZJRK8NaKrulSXLefCsu9wKShhokbH/dugJLwn9r0zoKd0mxpDbZMdGhLU1VVA+EomVrrqwE9RVChqZ9Smzgx73FCypKKslAunE0g7dZK0hHjSTp2kODeHjKTTVy3p5eDmjkuHjoZgztsHV28fbJzaSN05USt7e3vMzMywsrK6bjdpZWUlX3zxBb6+vgCMGzeO7777jgsXLmBjY0NAQACRkZFER0czYcIEnJycsLKyumYX7ObNmzl69Chnzpypzp4tXryYbt26ERMTQ1BQEBEREWg0Gl566SU0Gg3Dhg0jPj6e7du3M2LECDQaDS+//HKtx09ISECv1+Nfxyx+f39/cnNzyczMxNXVlcjISFauXAlAbGwsZWVlBAYGMmjQIDQaDVOnTkWj0WBhYVEj+wfwyiuvMG/ePEpKSggODmbNmjXX/H3WV4ODt3/84x88+uijvPvuu1hZyV9ztzUHTxjxHgyaAXu/Rrd7Ph3KMnhH9S3Pm/zEwqoRfK8bhrezXCfGYurmhuuLL+D8zNOkrFpG4bIVmJ1MwuzPgyT9OQGLHj1wfOgh7O66E5W5eVM3t9VrSFB3PSaWeu54sDPblifUGviZWVji2a0nnt16AobuuAPrj7Nz5Xa0VZnotRmYmuVSXpxHXnoaeelpJOy5PNbG3NoGO5f2uHh1oK2PN23aeeLUvj02jrc+qKutbl5rZWJuznOLfmyy925sVlZW1YEbQNu2bfH29sbGxqbGtoyMjHofMy4uDk9PzxrdngEBATg4OBAXF0dQUBDh4eEsWLAArVbLli1bGD58OG5ubmg0Gnr27MmpU6fq7JK9pL7dyBEREcyZM4e0tDQ0Gg0DBw5ErVYTHh7O/PnzAUM2LjQ0FPO//I5nzJjBtGnTOHv2LLNnz2bSpEmsWbPmpj9jDQ7ezp8/z3PPPSeBm7jMygkiXkEV+iwHfpmL6/H/0V7JYobpD7yg/g3TXYcg5Bmwb9/gQ8vkh/pRmZvT4cEp6B+YTNmRI+R+/z0Fa/+g7OhR0l59lYwPPsDh/nE4THgAs/btpKzILXat5bquxS/EjY49XOoV+BXnlbNndQYq0y6oTLsAF4sLz+5OSV4qmUmJZJw9Q+bZM2SnJFNeXERmcTyZSfHEbr18HDNLK0Mg184TmzbumJq3wb2TF+6dvTA1v3qd2NqCroYEY9eqm1eblh7oKYpSr67LluKvXYCKotS67a9DrG7WoEGDKCws5MCBA2zdupV3330XNzc33n//fXr16oWHhwedO3eu9bWdOnVCURTi4uJqrZIRFxeHo6MjLi4uAISFhWFmZkZ0dDTR0dHV3bFBQUFkZWWRmJiIRqOpddKms7Mzzs7OdOnSBX9/fzw9Pdm9ezchISE3df4NDt6ioqLYt28fPj4+N/XGohUys6bP+FdJy5lOwv6VeJ/4H6ZZsbD7v4bivz3uh7DnDd2u9bAiJplXVx1FpweVAu/d14MJQV5GPomWTVEULHv1wrJXL1xfeYW8lT+Su2I5ValpZH/9P7L/t4DC/l35rOMpDnvrUVRqZoXM4r7O9zV100Ud6hv41VXGpKLElA49etOhR2/AEPwsenUruqps9NpsdLoc9NpsbJ1KKchMp6K0hLRTJ0g7daKWtjjh4OaBfVs3KsttOHOkEkVlh0ptR8QjgXQb2L5Bwdi16ubVds4NDfTEjTMzM6ueZHCr+fv7k5KSQkpKSnX2LTY2lry8PAICAgBwcHCgZ8+ezJs3D1NTU/z8/HB1dWXChAmsWbOmzvFuAG3atGHYsGF8/vnnvPjiizXGvaWnp7N06VImTZpUnR2ztLRkwIABaDQatmzZwowZMwBD4BocHMyCBQtISUmpc7zbJZcC2PLy8mvuVx/1Ct4u1XUDuPvuu5kxYwaxsbH06NHjqgh75MiRN90o0bK5O9nBsGkw9FE4tdlQZiRpGxxeZnh0GWEI4rxCuGqK3kVp+aXVgRvI5IcbYeLkhPOTT9DmsWkUaTTkLl1K8c5d2O6J57U9kOoIG/vo+ajkLSkr0grUZ81WMAR5YILKpC2YtEV9cfud0wNp62NDXnoqqScSiV6yE11VDnpdHnpdHujLKcrNoSg3h3Nxx656/3WfqdnxvTPF+eYoKlsUlR2orNn8bSKW1iG4eLlhZe+ASq2u0Za6ZuP+NXhraKAnbo63tzd79uwhKSkJGxsbnJycbtl7Dx06lB49ejBx4kQ++eQTqqqqeOaZZwgPD6dfv37V+0VERPDZZ58xbpxhAoiTkxP+/v6sWLGC//73v9d8j3nz5hEaGkpUVBTvvPNOjVIh7dq1Y86cOTX2j4yM5OOPPwagT58+1dvDw8P58MMPqyc2XLJnzx5iYmIYOHAgjo6OnD59mjfffBNfX9+bzrpBPYO30aNHX7Xt0gyUKymK0mSRumiGFAU6DzU8zu03BHFxv8HJdYZH+/6GIK7rXYaSJFc4k1WM7i9f6jL54cYoajW2Q4ZgO2QIMXtWs/WzV4k4qscjFyZv1vGgRkdG7BvYTXkGy8BAGcTeQl0qY3K9yRHXCvJMTE1x9uxAWYkdJhY1P4B6XSmDJrTF1LyI5ONniN0eh16bh15XAPpiQEth9oVa27bqvV8M76OosLSzw9rRCRsHR8ys7KkqLQLFGkVlDYolKrUFJmalVFVYY3LFTMSGBHri5r300ktMnjyZgIAASktLOXPmzC17b0VRWL16NX/7298YNGgQKpWKESNG8Nlnn9XYLzw8nE8++aTG2LaIiAgOHz583fFunTt3Zt++fcyaNYvx48eTk5ODm5sbo0ePZtasWVcFq5GRkbz99tuMGDGiRrm08PBwZs2aRVRUVI1klpWVFatWrWLWrFkUFxfj7u7OiBEjeOONN64aF3cjGry2aUsma5s2A1mnYNdncOh70FYYtjl3gdDnoOd4MDFc1Gn5pYS9/2eNAE6tKGyfGdksgreWuvZkenE6UT9FYVqu5Y7jeoYe1OFzxf3WvHNnHB6YgP3IkahtbZuuobe5m7m+6lqz9UqxO1KvCvKu7H4syi1j8Ws7rwrwJs0xzHT96/N6vRYoYtAEdzRL9qHXFqDXFaDXF6PXFWNpU0lpQT56fcPGPZmYm2NpY4eFjQ2mFtZcOFMBWKCYuGJi3rNGmxpC1jYVzZFRF6ZvySR4a0YKL8Ce+RCzAMrzDdts3SH4aeg7FSzsWBGTzGurjqHV61ErCu/e173ZjHlrqcEbGJbSurQigwqF95yn0XtnBgW/r0VfZijyqlhaYnf3XThOeADLHk1XBf12dSuur+sFedcL8Op6vq7tOp2W0oICinJzKM7LoTg3l+LcHIrycsnPyKQwOwdtRQnlJUWUFRVdM9BTmXbE3G7MVW2qLwneRHNk1ODtueeeo1OnTjz33HM1ts+bN49Tp07xySefNLjBt4oEb81QWQEcWAS7PofCi2u+mdtBv0ch+GnSdPYkZZXg7WzVLDJul7Tk4A0MGbiUwhQ8bT2rx7ppCwrIX/0ruSuWU3Hqcp0wi27dDNm4u+5CZW1d5/Fk9mrjaS7X1/UCvLqer0/271r0Oh3lpSWUFRZSVlRIaVEhZYUFlBYVUZCZg4m5E72HD7vh7lIJ3kRzZNTgrV27dvz666/07VtzTcsDBw4wcuRIzp071/AW3yISvDVjVRVw9AfDyg1ZJw3b1GaG9VVDnwfnTk3bvr9oLjdXY9Dr9ZQeOEDu8hUUrluH/uJSLipra+xHjcRhwgQsunat3r9GJk9RyezVRtCar6/mQII30Rw15Fpp8Nqm2dnZ2NvbX7Xdzs6OrKyshh5OCAMTMwh8GJ7ZAw8sA88BhjFxBxbDvH6wfCKk7G3qVt4WFEXBqm9f2v37X3TaugXXGTMw7eCFrriY3O+XcWbUaJIeeJC8X34hLftsdeAGhjVVZ++aTXpxehOfhRBCtF4NDt46derEunXrrtr+xx9/SO03cfNUKvC7C6ZtgEfXQ5c7AT3Er4EFw+CbEXDiD2jkgo+idiaOjrSZ9ii+f/yB17ffYBsVBSYmlB46RNrMV8kdMZaHN1binn05ga/T60gpTGnCVgvR/NxGw8vFDWrINdLgIr1///vfefbZZ8nMzGTw4MGAYR2yjz76qFmPdxMtkFcwPLQcMk/AzrlweAUk7zI8apmhKoxHUamwDgnBOiSEqsxM8n5aRd4PP1CZmso9MXBPjJZYT4juqWKvvxpPW8/rH1SI24D6Yl27ioqKqxZBF+JKJSUlwNWrVtTmhmabfvHFF8yZM4fUVMMAc29vb956661mvyi9jHlr4QrSDDNU930D5QWGbTZuEPyUYYaqpcMta4qMSQK9Vkvx9u3EfvMJ1nvjUV38JtFamtPm3lE4jBuLRY8e16wbJxMdaifXl3HdynuBXq8nOTmZyspKPDw8UKka3OElWjm9Xk9JSQkZGRk4ODjg7u5+3dfcVKmQzMxMLC0tayxA25xJ8NZK1DZD1cwW+k2BAU+DfTujr4kqN9eaUhOPkvHTD1iu24nufGr1dvPOnXEYNxa7kSMxcXSs8RqZ6FA3ub6M61bfCyoqKjhz5kyjr+8pWhcHBwfc3NzqVShd6ryJlquqAo79CDvmQmacYZvKhDMed/F0YhjxOk+jrYkqN9fa6XU6SvbGkPfTTxRu2ID+4hp+iqkpNkOG4DB2LNahIVwoyyTqp6jqiQ4AKkXF+rHrJQOHXF/G1hT3Ap1OR0VFxS15L9HymJqaVnex10eDx7wJ0WyYmEHvh6DXg5Cw0VBm5Ox2Op77lXVmvxKt7cVX2nt4bRWyJuotoqhUWAcPwDp4ANo33yB/zRryf/yJsthYCteto3DdOkzc3Ske1p821loyHS7/hXlpooMEb6I1UqlUUipENBoJ3kTLpyjQZTh0Gc6RPZtJWfM+d6piiFQfJlJ9mGM6b4pinofIR0AtWYxbRW1nh9NDD+H00EOUxcWR9+NP5K9ZQ1VaGuaLV/MZcMxbIbqnwt4uClozmegghBD1ISMnRavi4hfK36peILLiIxZVDaNUb0Z3VRKdt78In/aGnfMMY+bELWXh74/bm2/QeesWPD78EKuQYFRAzyQ9z/+q4+vPtHyxryd2ceeuOV0+vTidvWl7pY6cEOK2JsGbaFXc7S15774enMOdWVVTGVgxj6NdngVrFyg4Bxteh4+7wYY3If98Uzf3tqMyN8f+nrvp8O23+G7ahMXjk9C1dcaqHBw37OPsw49wethwMud+RkVyco3XrkpYRdRPUUzbMI2on6JYlbCqic5CCCGaVr0mLMydO7feB/zrmqfNiUxYuH2k5ZfWXBO1sgyOrIBd8y4vv6Uyge7jIPRZcOvRoOPLgPLGo9fpKNm3j/zVqylctx5dcXH1c5Z9+mA/ahSl4X24c8PY22aCg1xfxiX3AtHS1St469ixY42fMzMzKSkpwcHBAYC8vDysrKxwdXUlMTHRKA1tDPKBFeh0kLABdn4GZ7df3u4TaQjifIcYxtBdh9xcjUNXWkrhps3kr15N8c6d1Stp6E1N2dWpii3dFY50VNCqDf+Pvon6hiC3oKZsslHI9WVcci8QLV29JiycOXOm+r+///57Pv/8cxYsWEDXi4tTnzhxgscff5wnn3zSOK0UorGoVNB1hOFx/oAhE3f8F0iMNjxcuxmCuO7jDLNZxS2lsrTE/t57sL/3HiovZFCw5jfyf/mF8oRThMZBaJyePCvY5a+wo7ua9jbtm7rJQghxyzW4zpuvry8//vgjgYGBNbbv37+fcePG1Qj0mhv5a0vUKvesYeWG/Yug8mKXna07DHiyzpUbJDNy6+j1esrj4ti38CNMNu3EvuTyc6bt22N3993Y33M35p07V29v6Ss3yPVlXHIvEC1dg0uFpKWlUVVVddV2rVbLhQsXGqVRQtxSjh1gxHsQ/jLsXwi750NhGmx6C7Z+CH0mwYCnDPuJW05RFCwCAhj4rwWk5Z0jbcs67LccoVKzg8pz58j+8kuyv/wS865dsbvnbnYGqHjz9FxZuUEI0Wo1OPN27733cv78ef73v//Rp08fwJB1e+KJJ2jXrh2//vqrURraGOSvLVEvl1Zu2PkZZMQatilq6DYaQp6Fdn0kM9IM6EpLKYqOJn/N7xRt2waVldXPxbWH7d1U7OmqUGSjrnViQ3POzsn1ZVxyLxAtXYODt8zMTCZPnsy6deuqv1SqqqqIiopi4cKFuLq6GqWhjUE+sKJB9Ho4vdkQxCVqLm/3voOq/k/x+8lK7rr7Hrm5NgPavDwKNmzg3M/LUR2Mq66BpFPguJdC57GT6THucUycnIDmv66qBG/GJfcC0dLd8NqmJ0+eJD4+HgA/Pz+6dOnSqA0zBvnAihuWdsQwueHYT6AzDBsoNHfHcsjLZPqO4Uyelo7O1rIEVxNLL07ngYXDCY7VEnZcR6cra/mqVFgN6A+RIUwsnke+1eWvvtrKjjRlZk6CN+OSe4Fo6WRheiEaIv887JmPfv+3KOWFAGTp7VhUNZxluqHMuC+MCUFeTdzI29uVWbW2+QpvlA6hQ8w5yo4dq95Hp8CxDgq7/RRiuijkWys1yo40dWZOgjfjknuBaOluKHg7d+4cv/76K8nJyVRUVNR47j//+U+jNa6xyQdWNJbKohwOffca7dI34KFkA1CmN+UX3R0MmToLF5/eTdvA21x6cTophSl42npWZ80qUlIoXL+e7N9/Qxt3snpfHXCyvULg+KfxuOs+chzVRP0Udc2CwMbOyknwZlxyLxAtXYODt82bNzNy5Eh8fHyIj4+ne/fuJCUlodfr6dOnD3/++aex2nrT5AMrGktlZSWfLvuD+bF67lLtYZrJH/RSXVGg2ncIhDxT76K/4tb6bdv/2Pv9J/SP19IpreZz2s4d+NEthT1dFc45U/3/71Jm7lpZucYK6iR4My65F4iWrsHBW//+/bnzzjuZPXs2tra2HD58GFdXVyZOnMiIESN4+umnjdXWmyYfWNFYKisr+f7ntcw+aIJOD6Cnn3KCx0zWEaWOQeHix8rFH4Kfhp7jwVTGwzUnl7Jz7YrNsdh5hMKNGynZt696VQeANEfY11nhQGc1Hz+3DsWk7qzcztSdjdbVKsGbccm9QLR0DQ7ebG1tOXToEL6+vjg6OrJ9+3a6devG4cOHGTVqFElJSUZq6s2TD6xoLJdursVte/Lm6ji0ej1qReHd+7ozwVcLe76Eg99BRZHhBVZtoN80CHoMbNs2beNFnapycij6809O/rIEi4MnMNVefk5lb095/27Ms97NIR+FUvPLGdUPwz/k5a0v19nV2tCMnARvxiX3AtHSNbhIr7W1dfU4N3d3d06fPk23bt0AyMrKatzWXWHkyJEcOnSIjIwMHB0dGTp0KB988AEeHh5Ge08hruf+vu2J9HcjKasEb2ery7NN73wfIl+FA98ZVm/IT4Gt/4Idn0CP+yH4GXDr3qRtF1czcXLCYdw4+o8bR9qF06RFr8M+5gTa7XvR5udjunEnLwJVKoj1Ujjoq3Cokxq9TlcjcAPQ6XWkFKY0akZOCCHgBjJvo0eP5u677+bxxx/npZdeYvXq1UyZMoVVq1bh6OjIpk2bjNLQjz/+mJCQENzd3Tl//jwvvfQSADt37qz3MeSvLdFYGpQZ0VZB/G+w63M4t/fy9o7hEDIdOg0zrLkqmi19VRWlhw5R+Gc0aet/xfx8zT9UVe08WOuWzgFfQ125SlMFlaJiyZ1LePiPh685+aE2knkzLrkXiJauwcFbYmIiRUVF9OzZk+LiYv7xj3+wc+dOOnfuzH/+8x86dLg1Swj9+uuvjB49mvLy8np/uckHVjSWG765ntsHu/4LsatBf7FPrk1nw7i4Xg+CmZVxGiwa1fnje8n8cz3W+05QdeBIjdUdKkwMWTn3yBHYhA1k6sl/XjVp5cqyJLWR4M245F4gWroWWectJyeHp59+mvPnz7N9+/Y69ysvL6e8vLz654KCAjw9PcnKypIPrLgplZWVbNy4kWHDht3YzTX/HKp9X6M6+B1KeQEAektHCgImEt9+PO6ePrjbWzRyq4Ux6EpKKNm9h5Jt2yjcugV9RmaN57Nt4UhHhSPeCsc7KBTYqvl91O+0tap77OOV11dOZQ7Jhcl42Xpd8zWi/goKCnB2dpbgTbRYNxS85eXl8eOPP3L69GlmzJiBk5MTBw4coG3btrRr184Y7QTglVdeYd68eZSUlBAcHMyaNWto06ZNnfu/9dZbzJ49+6rt33//PVZWkuEQTc9EW4pnzjZ8MzZgXZEBQKVezR+6/qS0HU779j5SaqQl0esxu5CB1cmTWCckYJmYiKqqqsYu+c626H0DKPXxocSnI9prBA/7yvexunQ1evQoKIyyHEU/837GPotWr6SkhIceekiCN9FiNTh4O3LkCEOHDsXe3p6kpCROnDiBj48Pb7zxBsnJySxevLjex5o5cyYffPDBNfeJi4vDz88PMEyIyMnJ4ezZs8yePRt7e3vWrFmDUsfNTTJvwlhuOvP2F2m5xcz59BMeVf/BAFV89fZy116oQ55C7z8K1GY3/T7i1tKVlVF24CCZWzZSGhODyalkw5q5VzD19sYyKAjLoH5YBgVh4uxMZWUlP67/kY8KPkJHzfFy18vaXSi5IJm665DMm2jpGhy8DR06lD59+vCvf/2rus6bj48PO3fu5KGHHmpQqZDMzEyys7OvuY+Pjw9mZlfftM6dO4enpyc7d+4kJCSkXu8n4xxEY2nsMUk7T2fx0Nd7AAhQkpiiXs8o9U7MlYtjqWzaGkqN9JsKNq43/X6iaWjz8ynZv5+SPXspjtlLeVz8VcGcmY8PFv36slXJ4b+uGvJt6j9erqmX9Wop5F4gWroGlwqJiYnhyy+/vGp7u3btSE9Pr+UVdXNxccHFxaWhTQBAd7GQ5pWZNSFaqo7O1qgU0OkhVu/Ny1VP8m/tQ/wZeQbbo4ugMA0078K2D6H7WBjwFHj0bupmiwZS29tjO3gwtoMHA7UHcxWJiVQkJtIb+Bo41wZOtFdI8FA41U5Fe6vayyOlF6dXB25gKFUye9dsQj1Ca53ZauwlvoQQxtPg4M3c3JyCgoKrtp88efKGA7Hr2bNnDzExMQwcOBBHR0dOnz7Nm2++ia+vb72zbkI0Z+72lrx3Xw9eW3WsuuDvS/eFYhv0AAx5yTA7dc98OBcDh5cZHp7BEPwU+N0L6gZ/lEUzcFUwl5dHyf79FO7azYU/N2Oelk77bD3ts/UMOawHdBQsG0VFjx5Y9uyJZe9eWPbsiYmzM8kFyXXWmvtrcFbfDJ0EeEI0Tw3uNn3sscfIzs7mhx9+wMnJiSNHjqBWqxk9ejSDBg3ik08+afRGHj16lOeff57Dhw9TXFyMu7s7I0aM4I033mjQBAlJlYvGYqxSDmn5pVcX/L3Suf2GIO74z6C72KVq1w6CHiO98wQSi83p6Gxd+2tFi3Hp+ooKCyNjj4bc/XuwOnkebewJ9CUlV+1v6uEB3buyQLeVE+5wti1UXKw199eacunF6XUu8XXlfvUJ8BoruLvVQaLcC0RL1+DgLT8/n3HjxrFv3z4KCwvx8PAgPT2dkJAQ1q5di7W1tbHaetPkAysaS5PX4SpMh5gFsO8bKDEUjC3Tm/KLNowluuE8MuZeJgR53fp2iUZR1/Wl12opP3Wa0sOHKD1yhLLDhyk/dfqqcXM6BVKdwLZbTzr3H4a5nz8W/n6YtGnD3rS9TNsw7ar3vHIsXX0CvGsFdw0JxppinJ7cC0RLd8N13rZv386RI0coKiqiT58+DB06tLHb1ujkAysaS5MHb9UNKSMvZjnn1n1Md1VS9eZ9uq743vMijn3GgonMUm1pGnJ9aYuKKDt6lNLDRyg9fJjiw4fQ5+TWuq+Jiwt07siP+n2ccYVkV4V0R9CbqGsEZtcL8K4V3DVkObD6ZgH/+pqbzdLJvUC0dDc8UGbgwIEMHDiwMdsihGgoUwti297DQxXO9FNOMMlkI3eq9tJPdQLWPgVbZ0HfKdB3Kti5N3VrhRGobWywDgnB+uL4X71eT1VmJuXx8ZTFxVMWH2eYCHH2LFWZmZCZyegrXl+lAq2HE1Ux75Lh64O5ry8eHg5YVCmUmVz+216lqPC09QSoc3zd4czDDZo00ZBxeiCzaYW45IaCt82bN7N582YyMjKqZ31e8s033zRKw4QQ9WOYqaqwT+/Hvko/XMhlosmf/M1uG+qiC7DlA9j2EfjfC/2fAK8QKfzbiimKgqmrK6aurtgMGlS9XVdcTNnJk9VBXeHxI1QlJmFSWobJuUwKz22EjZePs0hRyLDXk+aocMFJoXfv4VjviaO8QzGejm6oFNVVGTO9Xt+gYMzLzqvW41wKEq/U0Nm0QrRmDQ7eZs+ezdtvv02/fv1wd3evs0CuEOLW+OtM1RzFCfdRb6EOdIP432Dv15C8yzDJ4fjP0LY79H8cetwPZs13jKpoXCpra6wCA7EKDATAnYtZuvR0yk8nUnH6FOWnEylPPE3FqdNo8/Jomwdt8/RwRg/713JuwVrDwRSF75ztOWGVT6Y9ZNsphPYdSZckHe55Chm2OrRqw72hrmAMwM3ajVkhs67KpjVGlk6I1qzBY97c3d3517/+xSOPPGKsNhmNjHMQjaXZjHm7wjVnqqYdgZiv4chKqCoFoMrMjvJuE7AOfRxcujZBi0VdmsP1VZWTQ8Xp01QkJ1NxNpmK5GQqk5OpOHsWXXHxNV+rA/JsINdGwbV9Fzw79sDExeXqh7MzipkZ6cXppBSm4Gnrec2xbg0dH1cXuReIlq7BmbeKigpCQ0ON0RYhxE1wt7esu0SIe08Y+RkMnc2h3/6LY+xiOlRkYHLwazj4NXjfYVi9we9emeAgADBxcsLEyQmroJqrOej1erS5uVScPUtlSgqVqalUnk81/HvxoSovx6kInIr0kH6CvH0n6nwftYMDJi4utHVxQevgQLqDA2oHB9QO9hf/dUBtb4+TgwOze7zMW0c+QIv+mlk6IVq7Bgdvjz32GN9//z1vvvmmMdojhDCitEpL7jvUB72+N+GqI0xUb2aw6gDqpG2QtA2sXSDwYcMkB0fvpm6uaIYURakO7LjYBXslvV6PNieHytQ0qjIz635kZUFVFdq8PLR5eZQnJFz3vbsCy1QqqkJ64TL3PxK4idtWvYK3v//979X/rdPp+Oqrr9i0aRM9e/a8KqX/n//8p3FbKIRoNGeyitHpAVRodL3R6HrjTjY/9D+JZ+JKKEqH7R/D9k+g01Do9yh0Hi4rOIh6UxQFkzZtMGnT5pr76XU6tHl5hkAuwxDMafPzqoM5bX5+zX/z8g0FinU6HK2cJXATt7V6fSMfPHiwxs+9e/cG4NixYzW2y+QFIZq3K9dQvSRDccZkyP0wahac+MNQ+DcxGk5tNDzs2kGfydBnkpQbEY1GUakuZ/C61m/Mpa6iAm1eHtxQdVIhWo96BW/R0dHGbocQ4haobQ3Vd+/rfnmsXMBIwyP7NOxfCAeXQMF50LxrKDnid5chG9cxAlSqJjwTcTtSmZmhcnVt6mYI0eSkL0SI28yEIC8GdXG59hqqbXxh+P9B5OsQ96shG5e8C+J+MzwcOxomOPR+GKyv3T0mhBCicUnwJsRt6JozU69kagE9xxseF2Jh/7dweDnknoGN/4Q/34GA0YZsnFewFP8VQohbQPo9hBD10zYA7vo3/D0O7p0L7r1AWwFHf4BvR1DyST8K/vwPFGU2dUuFEKJVk+BNCNEw5jbQdzI8uRUejybR8z5K9OZY5Z/CbutsdB/5wYqHIWEj6LRN3VohhGh1pNtUCHHD0mz8GXpqHFb6u7hXvYsJ6mh6qxIvj42zawe9HzLUjpO6cUII0Sgk8yaEuGGX6sYVYcUy7RBGV7xDVPn7pPlNAUtHw0zVrf+GT3vBopFw9EeoLGvqZgshRIsmmTchxA2rrW7cKTrAnVPA+l8QvwYOfAeJGjizxfCwcICeE6DPI+DWo4laLoQQLZdk3oQQN+xS3Tj1xVmmNerGmZhD97Ew6Rd44QiEzwS79lCWB3u/hPkDqfj8Dk6v/ZT0C2lNeh5CCNGSSOZNCHFT6lU3zsELIl+F8JcNqzcc+A5t3BrMMo7gm3GE8j1vk+IeiWfkY9BpCKhNrz6GEEIIQII3IUQjqHfdOJUaOg0lzSWMew/+wijVdsapt+KvSsYzfSMs2wjWLtDjfuj1oKFbVWrHCSFEDRK8CSFuuTNZxWTp7VigvYsF2rvwV84yVr2VSTZ7MSvOhN2fGx6u3aD3g4ZgztaNtPxSzmQV09HZun7BohBCtEISvAkhbrm/TnSI03fgPe0k7n7ya9wzd8LhZRC/FjKOw4Y3YOM/SXUO5f3U3qzX9qNSMeO9+3owIciraU9ECCGagARvQohb7tJEh9dWHUOr11+e6OBoC45R0CUKSnPh+C+GQC5lDx6Z25lrup0CE0vWafvz289hDOo0HXdHm6Y+HSGEuKUkeBNCNInrTnSwdIR+U6HfVA4cjGHbT/9lrHob7ZUsxptsYTxbqPjqa+g5ztCt2q6PjI8TQtwWJHgTQjSZ+k50cPfpzqfa+/mkaix9lZOMVu/gbvUeHEszYc8XhodjR0MQ1+N+0sw8ZWycEKLVkjpvQohm71I3q0pRs0/vxyztY2y6ays89IMhYDO1gtwzsPVf8N8gsj8KJvqbN7n//R9YEZPc1M0XQohGJZk3IUSLUHs3q69hfFxFMZz4g7IDy1En/kl3VRLdVUm8zvcc+K0zBYUTset7P9i3b+rTEEKImybBmxCixaizm9XMGnqM44BVBM/EbeIu9V5GqnfSX4mnjyoBtr5leLTvD93GQMAo0nCSrlUhRIskwZsQotXo6GxNgWLL99ohfK8dggt53KWO4TXvOMzP74Fzew2P9a9yXteFTdoBbNAF8bf7IqXsiBCixZAxb0KIVuOva63mKI4EjP475o+vh3/Ew53/prxdMDq9Qj/VSf5p+h3bzZ/D/7dRFGz8ADJPNvEZCCHE9UnmTQjRqtRZgsTWDQY8wX7n+3jh63Xcpd7Dneq9BCkn6KlKhB3vGh7OXcH/XsPDvZeUHxFCNDsSvAkhWp1rlSDp6GxNluLIQu0IFmpH4Ew+w9X7+Wfn01gkb4esE7DtBGz7kHJrD6o6RWHd417wHggm5rf4TIQQ4moSvAkhbit/Xd0hV3Gg1+jnsQjygrJ8OLmBlJ0raJO2FaviVMwPfwuHvwUzW+g0BLreSXrbO0gsNpfJDkKIJiHBmxDitlNn16qFPWkd7iF8mRWm+omEqY4xVLWfoeqDuFbkQewvEPsLLnqFs/qufKfrTc/IcYwYPBQUhbT8UpnBKoQwOgnehBC3pbq6Vs9kFaPTQzlm/Knrw5+6PrxepePXMdZ4Z2/h/O6f8FOlMECJZ4AqHrYtR3vQnSTHUP6V6MV2bXdKFUN279IMVgnqhBCNSYI3IYS4Qkdna1QK6PSXt6kUNc5+IRzN6slDW/vTXskkQnWISNUhQlXHsSxKw7foJ740hQoTNft1Xdm2uieZDtP4M9eVV38+jk4PKoUaQZ0QQtwICd6EEOIKfx0Tp1YU3r2ve3XGTKXAOb0LS7TDWKIdhpVSybeR5cRu/YnBqoN0UGUQoo4lhFhYupwhejs+MunBdm0Pduv8eW3VMQZ1cZEMnBDihknwJoQQf1HXmLjaArtZ9/XBq4sLD0bbMLtqEt5KOoNURxikOkqEWTzO2gLGqHcwRr0DgGSdC+pfw6HHUMMMVgfPpjxVIUQLJMGbEELUoq4xcXUFdpeCuiS9Oyk6D7qNfolMX3te/PBr7lAdJkQVS08lES9VJpz+0fAAcPQ2BHHedxj+/cv6qzJeTgjxVxK8CSFEA9UW2NUV1I0ZM57XVgWgrdJjp5Tx2cByws1OQNJ2SD0IuUmGx8ElhgM5dkTtFUr7PBvW7GjLP9Zny3g5IUQNErwJIUQjaUhQB0B5ISTvhqRtF4O5Q5B7BlXuGfoCfc9+SaBpW/bpu3JQ14klPycxqNMk3B1tb+l5CSGaFwnehBDCyOpc8cHcFjoPMzwAygogZQ/a0xoyDq7FtewM3qoLeHOBceqtAGjn/R+0C4T2/S4+gsDOo/qQ0s0qROsnwZsQQjQXFnbQeRg67wg2FPTno4MV9FFOEKg6RaByit6q09hpSyB5p+Fxia0HtO/HIX0n3jtiw1GdN2WKhXSzCtFKSfAmhBDNkIM5zBzVjzdXW6OpCjSULBkTwISO5XBuH5yLgfP74MJxKEyFuF/pDawwA51eIUnflrhfO1CQNxi7Dn3ArQfYuoGiNPWpCSFukgRvQgjRTN3ftz2R/m5Xj5dz6QqBEw3/XVEMqYdIOqIhPiaa3qpTuCm5+Cjp+JAOO/bAjosHtHI2BHFu3cGtJ7gGQJtOYGrRJOcnhLgxLS54Ky8vZ8CAARw+fJiDBw/Su3fvpm6SEEIYTZ3j5S4xswbvMMwd+/DMrt7o9NCGfPxVyXRXneX57mVYZsdC1kkoyYLEaMPjIr2iQnH0Bhc/cO5iCAxduhr+21wmRgjRHLW44O3ll1/Gw8ODw4cPN3VThBCi2biygHC23p5d+p7cO+ohLC+NeasshYxYSD9KwpGd5J85SBflHHaUQE6i4XFibY1jam3bUWjdAVOXTli7dQInH8PDsSOYWTXBWQohoIUFb3/88QcbNmzgp59+4o8//mjq5gghRLNyzbIkppbQri9pNgFE/eiCTj8K0ONCHn6qVOZFWWFfdAYyTxgexRmoC8/jUHge0nfC0b+8ma27IYhz8DKsEmHvaSgw7OBl+NdUZroKYSwtJni7cOECjz/+OL/88gtWVvX7i6+8vJzy8vLqnwsKCgCorKyksrLSKO0Ut4dL149cR8IYbub6crYywdnLrs7Xn0ovQKe/9JNCJo5k6hw56t6PAR2dAEjLL2PUR2vx4TzeygU6qNLpqFxguHsJZgVJKGX5UJhmeFw56/UKemsXsHFDb+0KNm3R27Q1/GvbtsbPmNz68XbyuRUtXYsI3vR6PVOmTOGpp56iX79+JCUl1et17733HrNnz75q+4YNG+odAApxLRs3bmzqJohWzBjXV145KKjRc3nWqYKe04d2kx1n+DkhXyFXb8N+urJf3xV0hu3P2mjp3E6PaVURiWkZnEzNpJ2STTslk37WmXgo2VhVZGGiK0cpzoTiTK43t7VCbUW5qQNlJg5UmNgYHmobKk2sqVBf+fOl/7YGRXVTv4OSkpKber0QTU3R6/X66+9mHDNnzuSDDz645j5xcXFs2LCBH374gS1btqBWq0lKSqJjx47XnbBQW+bN09OTrKws7OzsGus0xG2osrKSjRs3MmzYMExNTZu6OaKVMfb1tXL/Od5YHVu97NY7owK4v+/lNVXT8suI+GjrFRk6w36afwzC3d7ims9vT8jk37/G4E4WbVV5PBloRVCbCii6gFJ0oea/2nJuhK7TMLQTlt3o6VNQUICzszP5+flyLxAtUpNm3v7xj38wZcqUa+7j4+PDn3/+ya5duzA3N6/xXL9+/Zg4cSKLFi2q9bXm5uZXvQbA1NRUbriiUci1JIzJWNfXQ8Eday9BcpGXs2n15AetXm+oMXdfd7ycDbNPz+Xn1wjcAHR6OHK+kDd+jUOntyEXG2K1sHW/wvaZkVfPmNXroSyPP3Yd5PtNe2lDPo6qIsb6W9HdQQuluVCaY/i35OK/5YahL7kVKipKqm54BQn5zIqWrkmDNxcXF1xcXK6739y5c3nnnXeqf05NTSUqKooVK1YwYMAAYzZRCCFapeuVILnW5IeOztaoFGoEcGpFgb9sA9Dq9SRllVz9XopCWoUF0zeWotP3MGzTweJjdQR7wA97TvPvX/bAST3Z7/8pK0iI21aLGPPm5VXzw2ljYwOAr68v7du3r+0lQgghblJdAd6VZUmuzMz17eBYa1Dn7Vz7GOMzWcX1DvbS8kuZ+Us8Or29YYMeXlt1jEFdXGQNV3HbaRHBmxBCiOalrsxcbUFdXcFVXRm82oK9hgR6QrR2LTJ48/b2pgnnWQghhKD2zNw1a83V8vr6BnsNCfSEaO1aZPAmhBCi+brukl5XqG+w15BAT4jWToI3IYQQTaq+wV5DsnpCtGYSvAkhhGgxGpLVE6K1urky1UIIIYQQ4paS4E0IIYQQogWR4E0IIYQQogWR4E0IIYQQogW5rSYsXKoNV1BQ0MQtES1dZWUlJSUlFBQUyDqJotHJ9WVcl+4BUi9UtFS3VfBWWFgIgKenZxO3RAghRFMrLCzE3t6+qZshRIMp+tvoTw+dTkdqaiq2trYoitLUzalTUFAQMTExLe69bvRYDX1dQ/a/3r43+nxBQQGenp6kpKRgZ2dXr7Y0F3J9Nc7+9dlPrq/m+V56vZ7CwkI8PDxQqWT0kGh5bqvMm0qlahEL2avV6lv2hd2Y73Wjx2ro6xqy//X2vdnn7ezsWtzNVa6vxtm/PvvJ9dV830sybqIlkz85mqHp06e3yPe60WM19HUN2f96+97s8y2RXF+Ns3999pPrq/W8lxDNyW3VbSpEYykoKMDe3p78/PwWlxkRzZ9cX0KIa5HMmxA3wNzcnFmzZmFubt7UTRGtkFxfQohrkcybEEIIIUQLIpk3IYQQQogWRII3IYQQQogWRII3IYQQQogWRII3IYQQQogWRII3IYQQQogWRII3IYwsJSWFiIgIAgIC6NmzJytXrmzqJolWZMyYMTg6OjJu3LimbooQ4haRUiFCGFlaWhoXLlygd+/epKf/fzt3ExJV/4Zx/JqmhlLDCsXRsgaCJAua0ByMhAzDNkGR0UJSQ1wEWVEWQoUuwoI2FprYolHcTBHVosyF0lSGvRmRvZBlqUS+IBaVlsWMzyKe+T/+J6VMHU9+PzCLuT33ObfDb3HxO2emU3FxcWpublZwcHCgR8NfwO1269OnT6qsrNSFCxcCPQ6ACcDOGzDOIiMjZbfbJUlWq1VhYWHq7e0N7FD4a6xdu1azZ88O9BgAJhDhDVPezZs3tXHjRkVFRclkMuny5ct+x5SWlspms2nmzJlyOBy6d+/eqK7V2Ngoj8ej6OjoP5waRjCRawvA1EF4w5TX19enFStWqLS09Kd/P3funPbt26eCggI9fPhQK1asUGpqqrq7u33H2O12LV++3O/17t073zG9vb3KyMjQmTNnxv1/wuQwUWsLwNTCM2/Af5hMJl26dEmbNm3y1RwOh1atWqWSkhJJktfrVXR0tHJzc5Wfn/9L5x0YGND69euVk5Oj7du3j8fomOTGa21JP557Kykp4Zk3YIpg5w0Ywbdv39TY2KiUlBRfbdq0aUpJSVFDQ8MvnWNwcFBZWVlat24dwQ0+Y7G2AExNhDdgBD09PfJ4PIqIiBhSj4iIUGdn5y+d4/bt2zp37pwuX74su90uu92upqam8RgXBjIWa0uSUlJStHXrVlVXV2vBggUEP2AKmB7oAYC/3Zo1a+T1egM9Bv5StbW1gR4BwARj5w0YQVhYmMxms7q6uobUu7q6ZLVaAzQV/gasLQCjRXgDRmCxWBQXF6e6ujpfzev1qq6uTomJiQGcDEbH2gIwWtw2xZT3+fNnvXr1yvf+zZs3evTokebNm6eFCxdq3759yszMVHx8vBISElRcXKy+vj7t2LEjgFPDCFhbAMYDPxWCKc/tdis5OdmvnpmZqYqKCklSSUmJTpw4oc7OTtntdp06dUoOh2OCJ4XRsLYAjAfCGwAAgIHwzBsAAICBEN4AAAAMhPAGAABgIIQ3AAAAAyG8AQAAGAjhDQAAwEAIbwAAAAZCeAMAADAQwhsAAICBEN6AScztdstkMunDhw8BuX5dXZ2WLl0qj8cz7DGFhYWy2+2+9/n5+crNzZ2A6QBgaiK8AZPE2rVrtXfv3iG11atXq6OjQ6GhoQGZ6eDBgzp8+LDMZvMv9+Tl5amyslKvX78ex8kAYOoivAGTmMVikdVqlclkmvBr19fXq6WlRVu2bPmtvrCwMKWmpqqsrGycJgOAqY3wBkwCWVlZunHjhk6ePCmTySSTyaTW1la/26YVFRWaM2eOrly5opiYGAUFBSktLU39/f2qrKyUzWbT3LlztXv37iG3OgcGBpSXl6f58+crODhYDodDbrd7xJlcLpfWr1+vmTNnDqkfP35cERERmj17trKzs/X161e/3o0bN8rlcv3x5wIA8Ed4AyaBkydPKjExUTk5Oero6FBHR4eio6N/emx/f79OnToll8ulmpoaud1ubd68WdXV1aqurlZVVZXKy8t14cIFX8+uXbvU0NAgl8ulx48fa+vWrdqwYYNevnw57Ey3bt1SfHz8kNr58+dVWFiooqIiPXjwQJGRkTp9+rRfb0JCgt6+favW1tbRfSAAgGFND/QAAKTQ0FBZLBYFBQXJarWOeOz3799VVlamxYsXS5LS0tJUVVWlrq4uhYSEKDY2VsnJybp+/bq2bdum9vZ2OZ1Otbe3KyoqStKP59JqamrkdDpVVFT00+u0tbX5jv9XcXGxsrOzlZ2dLUk6evSoamtr/Xbf/u1ra2uTzWb77c8DADA8dt4AgwkKCvIFN0mKiIiQzWZTSEjIkFp3d7ckqampSR6PR0uWLFFISIjvdePGDbW0tAx7nS9fvvjdMn3+/LkcDseQWmJiol/vrFmzJP3YJQQAjC123gCDmTFjxpD3JpPppzWv1ytJ+vz5s8xmsxobG/2+NfrfwPf/wsLC9P79+1HN2NvbK0kKDw8fVT8AYHiEN2CSsFgsI/6e2mitXLlSHo9H3d3dSkpK+q2+Z8+eDaktXbpUd+/eVUZGhq92584dv94nT55oxowZWrZs2egHBwD8FLdNgUnCZrPp7t27am1tVU9Pj2/n7E8tWbJE6enpysjI0MWLF/XmzRvdu3dPx44d09WrV4ftS01NVX19/ZDanj17dPbsWTmdTjU3N6ugoEBPnz71671165aSkpJ8t08BAGOH8AZMEnl5eTKbzYqNjVV4eLja29vH7NxOp1MZGRnav3+/YmJitGnTJt2/f18LFy4ctic9PV1Pnz7VixcvfLVt27bpyJEjOnjwoOLi4tTW1qadO3f69bpcLuXk5IzZ/ACA/zENDg4OBnoIAJPTgQMH9PHjR5WXl/9yz7Vr17R//349fvxY06fzZAYAjDV23gAM69ChQ1q0aNFv3cLt6+uT0+kkuAHAOGHnDQAAwEDYeQMAADAQwhsAAICBEN4AAAAMhPAGAABgIIQ3AAAAAyG8AQAAGAjhDQAAwEAIbwAAAAZCeAMAADCQfwBok/WrL0UaywAAAABJRU5ErkJggg==",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# plot the fitted curves\n",
"hm1_1 = ml.head(r1, 0, t1)\n",
"hm2_1 = ml.head(r2, 0, t2)\n",
"hm3_1 = ml.head(r3, 0, t3)\n",
"plt.semilogx(t1, h1, \".\", label=\"OW1\")\n",
"plt.semilogx(t1, hm1_1[0], label=\"timflow OW1\")\n",
"plt.semilogx(t2, h2, \".\", label=\"OW2\")\n",
"plt.semilogx(t2, hm2_1[0], label=\"timflow OW2\")\n",
"plt.semilogx(t3, h3, \".\", label=\"OW3\")\n",
"plt.semilogx(t3, hm3_1[0], label=\"timflow OW3\")\n",
"plt.xlabel(\"time (d)\")\n",
"plt.ylabel(\"head change (m)\")\n",
"plt.legend(bbox_to_anchor=(1.05, 1))\n",
"plt.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"### Comparison of Results "
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"The aquifer response obtained with `timflow` is compared to the results based on Hantush’s analytical solution (Hantush, 1955), implemented in the software AQTESOLV (Duffield, 2007). The results are almost identical."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
" \n",
" \n",
" | | \n",
" k [m/d] | \n",
" Ss [1/m] | \n",
" c [d] | \n",
" RMSE [m] | \n",
"
\n",
" \n",
" \n",
" \n",
" | timflow | \n",
" 224.63 | \n",
" 2.13e-04 | \n",
" 43.88 | \n",
" 0.060 | \n",
"
\n",
" \n",
" | AQTESOLV | \n",
" 224.72 | \n",
" 2.12e-04 | \n",
" 44.00 | \n",
" - | \n",
"
\n",
" \n",
"
\n"
],
"text/plain": [
""
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = pd.DataFrame(\n",
" columns=[\"k [m/d]\", \"Ss [1/m]\", \"c [d]\", \"RMSE [m]\"],\n",
" index=[\"timflow\", \"AQTESOLV\"],\n",
")\n",
"\n",
"t.loc[\"AQTESOLV\"] = [224.723, 2.124e-4, 43.998, \"-\"]\n",
"t.loc[\"timflow\"] = np.append(cal.parameters[\"optimal\"].values, cal.rmse())\n",
"\n",
"t_formatted = t.style.format(\n",
" {\n",
" \"k [m/d]\": \"{:.2f}\",\n",
" \"Ss [1/m]\": \"{:.2e}\",\n",
" \"c [d]\": \"{:.2f}\",\n",
" \"RMSE [m]\": lambda x: \"-\" if x == \"-\" else f\"{float(x):.3f}\",\n",
" }\n",
")\n",
"t_formatted"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n",
"\n",
"* Duffield, G.M. (2007), AQTESOLV for Windows Version 4.5 User's Guide, HydroSOLVE, Inc., Reston, VA.\n",
"* Newville, M.,Stensitzki, T., Allen, D.B. and Ingargiola, A. (2014), LMFIT: Non Linear Least-Squares Minimization and Curve Fitting for Python, https://dx.doi.org/10.5281/zenodo.11813, https://lmfit.github.io/lmfit-py/intro.html (last access: August,2021)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.13.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}