10_K.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 10.\n",
    "### Sensitivity of hydraulic head at a point to <font color='blue'>heterogeneous hydraulic conductivity</font> under steady state flow conditions\n",
    "\n",
    "## 0. Forward model\n",
    "<br>\n",
    "Governing equation:\n",
    "\n",
    "$\\begin{align}\n",
    "\\quad\\quad K\\ b\\ \\frac{d^2h}{dx^2}+R=0\\cr\n",
    "\\cr\\end{align}$\n",
    "\n",
    "Boundary conditions:\n",
    "\n",
    "$\\begin{align}\n",
    "\\quad\\quad -K\\ b\\ \\frac{dh(x)}{dx}=0\\ ,\\quad\\quad\\;\\;\\, x=0=\\Gamma_2\\cr\n",
    "\\end{align}$\n",
    "\n",
    "$\\begin{align}\n",
    "\\quad\\quad h(x)=h_{\\Gamma_1}\\ ,\\quad\\quad\\quad\\quad\\quad x=L=\\Gamma_1\\cr\n",
    "\\cr\\end{align}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import HTML, display\n",
    "def set_background(color):\n",
    "    script = (\n",
    "        \"var cell = this.closest('.code_cell');\"\n",
    "        \"var editor = cell.querySelector('.input_area');\"\n",
    "        \"editor.style.background='{}';\"\n",
    "        \"this.parentNode.removeChild(this)\").format(color)\n",
    "    display(HTML('<img src onerror=\"{}\">'.format(script)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "from warnings import filterwarnings\n",
    "filterwarnings(\"ignore\", category=DeprecationWarning)\n",
    "\n",
    "import numpy as np\n",
    "R, b, L, BC1h, ocol = 1e-1/1000., 10., 10000., 0., 5000\n",
    "X = np.arange(L)\n",
    "K1, K2 = 1., 10.\n",
    "K = np.hstack([K1*np.ones([1,round(2*L/3)]), K2*np.ones([1,round(L/3)])])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Direct sensitivity\n",
    "\n",
    "$\\begin{align}\n",
    "\\quad\\quad\\rm{Not\\ available}\n",
    "\\cr\\end{align}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Perturbation sensitivity\n",
    "$\\begin{align}\\cr\n",
    "\\frac{\\partial h(x')}{\\partial K}\\approx\\frac{h(x,K+\\Delta K)-h(x,K)}{\\Delta K}\\;\n",
    "\\cr\\cr\\end{align}$\n",
    "\n",
    "### 2a. Analytical\n",
    "$\\begin{align}\n",
    "\\quad\\quad\\rm{Not\\ available}\n",
    "\\cr\\end{align}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2b. Semi-analytical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAACdCAYAAACpbTGvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlyUlEQVR4nO3de1xVVfr48c+DqAdUUEGNEEXRQrmrKNaQGqZmajl2mbK8dHGaaSadpuz2bbLLfPNbTpaNNk2WWlmOZs10cdIynXJSE8W8ZF4IMswLl8ALoAjP7w+O5wcCCgocOD7v1+u8PGftvdd+2HJ49lp77b1EVTHGGGOM5/BydwDGGGOMqV2W3I0xxhgPY8ndGGOM8TCW3I0xxhgPY8ndGGOM8TCW3I0xxhgP4+3uAGpLYGCghoaGujsMY4wxpt5s3LgxS1XbnV7uMck9NDSU5ORkd4dhjDHG1BsR+aGycuuWN8YYYzyMJXdjjDHGw3hMt3xt++STT5g5cyZeXl7lXi+99BKdOnVi2bJlLFiwgCZNmpRb/pe//IWAgAA+/fRTVqxYQfPmzWnWrBnNmjWjefPm3H333fj4+JCSksLOnTtdy3x8fGjRogV9+/bFy8uLI0eOUFJSQosWLfD2tv8mY4wx1WdZowpFRUUcPnyYkpKScq+TJ08CkJmZyTfffFNh+fHjxwHYuHEjs2fP5sSJExQXF7vqvf322/Hx8eHtt99mxowZle7Xy8uLBx98kJdffhmApk2b0qJFC9q2bUtqaioATz/9NF988QV+fn74+fnh7+9PUFAQU6dOBWDdunUcPnwYf39//P39Xeu0aNGiTo+bMcYY9xNPmTimT58+2lAH1BUXF3PixAlOnDiBn58fIkJWVhZZWVkcP36c48ePU1hYSH5+PsOGDQPg888/Z/PmzRw7doz8/HyOHTsGwKxZswB44okn+Pe//82RI0fIy8vj8OHDtG3blvT0dACGDRvG8uXLy8URHh7Ojh07ABg3bhy7d+8mICDA9erZsyd33nknAMnJyXh7e7uW+fj4ICL1cbiMMcZUk4hsVNU+FcotuXuOkpISvLxKh1GkpqZy4MABV+LPy8vD19eX2267DYCpU6eSkpJCdnY22dnZ5OTkEB8fz+effw7AJZdcwu7du111OxwOxowZw1tvvQXAgw8+SJMmTejQoYPrFRYWRqdOner5pzbGmAuXJXdzVsXFxTRp0gSAtWvXsn//flfyz87Opnv37kyaNAmAbt26kZ6eXuGSw2uvvYaq0rVrV9q0acNFF11EcHAwwcHBJCUlkZiYiKqSk5ND27ZtrTfAGGPOQ1XJ3a65G5dTiR2gf//+Z1x3z549lJSUkJOTw8GDBzl48CABAQEAnDhxgkGDBnHw4EEOHDjApk2bOHToEKpKYmIimZmZdOjQAYfDwcUXX0zHjh0JDg5mwoQJDBkyhIKCArZs2UJoaCjt27e3EwBjjKkhS+7mnHl5eREYGEhgYCARERGu8ubNm/P666+XW/fEiRMUFRUB0KxZM2bOnMm+ffvIyMhg3759rF+/nqFDhwKwY8cOEhISgNLLAZ06dSI0NJT/+Z//ITExkezsbL799ltCQ0O5+OKLy52UGGOMseRu6smpW/4AWrduzZQpU6pcNywsjA8//JAffviB9PR017+nLiF9+eWXjB49GgBvb29CQkLo3Lkzs2bNIioqip9++ol9+/YRFhZG27Zt6/xnM8aYhsaSu2lw/P39GTFiRJXLExMT+eSTT8ol/x9++AGHwwHA0qVLuffeewFo06YNYWFhdOvWjRdffJH27dtz8OBBVJUOHTpYl78xxiNZcjeNTkBAgKsLvzK//OUv6dSpE3v27CE1NZU9e/awfv16fH19AXj++ed59tlnadGihSvxh4eH88QTT+Dt7c3JkyftwUHGmEbNRsubC05KSgpfffUVe/bsYc+ePezevZvc3FwOHDgAwC233MKqVasIDw93vSIjI0lKSnJz5MYYU56NljfGKS4ujri4uHJlZW/pGzp0KM2bN+e7775j0aJF5ObmEhMTw+bNmwG48847OXr0KD169CAqKorIyEjCwsJsYJ8xpsGw5G4M5W8DHD9+POPHjwdAVcnMzCQ7O9u1/Pjx42zYsIHFixe7BvmNHDmSDz74AIC//e1vdOrUicjISEJCQuy6vjGm3lm3vDHnKD8/nx07drBt2zYCAgIYMWIEBQUFtGzZkpKSEgBatWpFREQEv/nNbxg3bhyqSnZ2NoGBgW6O3hjjCaxb3pha5uvrS+/evendu7erzMfHh+zsbLZv3862bdtcr1P3+KelpREWFkbHjh2JjY0lNjaWuLg4EhMTadeunbt+FGOMh7GWuzH16NChQ7z55pukpKSwefNmvvvuO4qLi1myZAnXX389W7ZsYe7cua7EHxERQfPmzd0dtjGmgbKWuzENQPv27fnjH//o+lxQUMD27dsJCwsDYOfOnbz++uuuWQC9vb3p2bMnS5cupVu3buTm5tK8eXN8fHzcEr8xpnHwqs+diUgTEUkRkY+cn9uKyKcistv5b5sy6z4sIntEZKeIVH1TszGNmI+PD3369KFNm9Jf/RtuuIHDhw+za9cuFi9ezNSpUwkODqZDhw4APPfcc/j5+REXF8ddd93FK6+8wqZNm1zX+I0xBuq5W15E7gP6AH6qOkJEngVyVHW6iDwEtFHVB0WkJ/AO0Be4GPgMuERVi6uq27rlzYXgq6++YtmyZWzYsIHk5GRycnLw9/cnJycHLy8vFi5cyIkTJ4iPjyc8PNwexmOMh3N7t7yIdASuAf4M3OcsvhYY6Hy/AFgNPOgsX6Sqx4E0EdlDaaJfW1/xGtMQXXbZZVx22WVA6W16aWlp/PDDD3h5lXbCzZ49m7VrS78mvr6+9OrVi2uvvZb777/fbTEbY+pffZ7WvwBMBVqVKeugqvsBVHW/iLR3lgcD68qsl+EsM8Y4iQhdu3ala9eurrI1a9awa9cukpOT2bBhA19//TXbt28HSk8GLrvsMkJDQ0lISKB///7Exsa6JvQxxniOeknuIjICOKSqG0VkYHU2qaSswvUDEZkETALo1KnT+YRojEfw8vJyPTL31ltvBXA9aKegoICQkBDWrFnDokWLgNLpef/3f/+X++67j5MnT3Lw4EGCg+082pjGrr5a7pcDo0RkOOAA/ETkLeCgiAQ5W+1BwCHn+hlASJntOwI/nV6pqv4d+DuUXnOvyx/AmMbq1BPyfH19Wbx4MQAZGRmsW7eOtWvXEh0dDcA333xDnz59CAkJ4bLLLiMxMZHExEQiIyNd3f7GmMah3u9zd7bc73cOqHsOyC4zoK6tqk4VkQjgbf7/gLqVQHcbUGdM3dm/fz+LFy9m7dq1rFmzhn379gGwYsUKrrrqKtLS0ti/fz+9e/e2e++NaSDcPqCuCtOBxSJyB7AXuAFAVbeLyGLgW+AkcM+ZErsx5vwFBQUxefJkJk+ejKqSnp7OmjVr6N+/PwDz58/nySefxOFw0LdvX1fLPikpyUblG9PA2BPqjDHVkpWVxZdfful6paSk4HA4yM3Nxdvbm/fffx+AAQMG0LZtWzdHa8yFoaqWuyV3Y8w5OXr0KDt37nQ9W79///6sW7cOESEmJoYrr7yS4cOHk5SU5OZIjfFcVSV3GyVjjDknLVu2LDdpzn/+8x++/PJLnnjiCVq3bs3s2bN55ZVXXMufeeYZli9f7nq0rjGm7ljL3RhTJwoLC/n5558JCgri0KFDdOzYkaKiIry9venbty+DBg3i5ptvJiIiwt2hGtNoWcvdGFOvHA4HQUFBQOmEOT///DPLly/n/vvvp7i4mOnTp5OSkgLA999/z/Tp09m4caM9J9+YWmAtd2OMWxw+fBhvb298fX154403GD9+PACBgYEkJSVx1VVXceONN9KqVauz1GTMhcsG1BljGrQDBw7w2WefsWLFCj799FMOHTpEVlYWbdq04bPPPqOgoICBAwdasjemDEvuxphGQ1X5/vvvXfPcX3311XzyySd4e3uTkJDAkCFDGDZsGPHx8W6O1Bj3smvuxphGQ0RciR3g/fffZ+XKldx///0UFBTw+OOPM3XqVNfy5cuXk5WV5Y5QjWmQrOVujGl0srKyyMrKIjw8nNzcXAIDAykpKaFv374MHz6c4cOH06tXL3smvvF41nI3xniMwMBAwsPDAfDz82Pt2rU8/vjjqCrTpk0jPj6e2bNnA5Cfn09OTo47wzWm3lnL3RjjUTIzM1m+fDmJiYl07tyZd955h1tvvZWEhASGDx/ONddcQ0xMjGu2PGMaMxtQZ4y5IO3cuZOFCxeybNkyNm7cCEBISAgpKSkEBAS4OTpjzk+tJXcRaQEUNrRZ2iy5G2PO5sCBAyxbtozk5GTmzJkDwN133012djajRo1i+PDhlvBNo3LO19xFxEtEbhGRj0XkEPAdsF9EtovIcyLSvS4CNsaY2nbRRRdx++23uxI7lF6z/+9//8u4ceNo3749AwYMYMGCBW6M0pjzV50BdauAMOBh4CJVDVHV9kAisA6YLiK31mGMxhhTZ5599lkyMjL4+uuveeSRR8jNzWXr1q0AnDx5kkcffZQ1a9ZQXNygOiuNOaOzdsuLSFNVLTrfdeqadcsbY2rLyZMn8fb2ZvPmzfTt25eioiLatWvHtddey+jRo0lKSqJ58+buDtOYc++Wr07SdndiN8aY2uTt7Q1AbGwsmZmZLFq0iMGDB/OPf/yDa665hq+++gooHZl/5MgRd4ZqTKWqfZ+7iPQRkfdFZJOIbBGRrSKypS6DM8YYd/P39+emm27i7bffJjMzk2XLlpGYmAjA9OnTadeuHSNHjmTevHn2lDzTYHjXYN2FwAPAVsDmZDTGXHCaN2/O1Vdf7fp80003UVJSwnvvvcdHH32El5cXEydOZO7cuW6M0piaPaEuU1U/UNU0Vf3h1KvOIjPGmAaub9++zJw5k/T0dJKTk7npppvo0KGDu8Mypvr3uYtIEnAzsBI4fqpcVd+rm9BqxgbUGWOMudBUNaCuJt3yE4FwoCn/v1tegQaR3I0xxt2Kioo4efIkPj4+7g7FXOBqktxjVDWqziIxxphGbsqUKSxevJjMzEx3h2IucDW55r5ORHrWWSTGGOMBPGW+DtO41aTl/gtgvIikUXrNXQBV1eg6icwYYxoZEbHkbhqEmiT3YXUWhTHGeACbRtY0FGdN7iIiWqrK295OrVO7oRljTONjfwpNQ1CdlvsqEVkK/EtV954qFJFmOLvqKZ1cZn6dRGiMMY3E0KFDueiii9wdhjHVmjjGAdwOjAW6ALmAA2gCrABmq+rmOo2yGuw+d2OMMReac77PXVULgTnAHBFpCgQCBaqaW+tRGmNMI3bkyBEKCwtp166du0MxF7ia3AqHqhap6v6aJnYRCRGRVSKyQ0S2i8hkZ3lbEflURHY7/21TZpuHRWSPiOwUkaE12Z8xxrjDY489Rvfu3d0dhjE1S+7n4STwR1XtASQA9zjvmX8IWKmq3Sl9rO1DAM5lvwIiKB2lP0dEmtRTrMYYc85sQJ1pCOoluTtb+5uc748AO4Bg4FpggXO1BcB1zvfXAotU9biqpgF7gL71Easxxpwru8/dNBTVvs9dRK6kdFBdLrAN2AJsU9XjZ9quknpCgThgPdBBVfdD6QmAiLR3rhYMrCuzWYazzBhjGiy7z900FDVpub8FfERp0u0K/AnYXpOdiUhLYCkwRVUPn2nVSsoqnA6LyCQRSRaRZHuWszGmIbCWu2kIavKEuj2q+r7z/ZKa7sg50n4psLDMNLEHRSTI2WoPAg45yzOAkDKbdwR+Or1OVf078HcovRWupjEZY0xtGjFiBB07dnR3GMZU6z73N4BNQBBwSFX/UuOdlPZVLQByVHVKmfLngGxVnS4iDwFtVXWqiEQAb1N6nf1iSgfbdVfV4qr2Yfe5G2OMudCcz3zuC4AY4CJgqIjcC3zjfG1R1eq04i8HbgO2ishmZ9kjwHRgsYjcAewFbgBQ1e0ishj4ltKR9vecKbEbY0xDkJ2dzeHDh+nSpYu7QzEXuLO23CtsIOIN9KQ04Uer6gN1EVhNWcvdGONuDzzwAHPmzOHYsWPuDsVcIM6n5V6Oqp6kdKT8ltoIzBhjPIkNqDMNQX09xMYYYzye3QpnGgpL7sYYU4us5W4aAkvuxhhTS6zlbhqKGl9zN8YYU7nrrruOrl27ujsMYyy5G2NMbenfvz/9+/d3dxjGWLe8McbUlv3797Nli91IZNzPkrsxxtSSl156iT59KtxybEy9s+RujDG1yEbLm4bAkrsxxtQSGy1vGgpL7sYYU4us5W4aAkvuxhhTS6zlbhoKuxXOGGNqyZgxY7j00kvdHYYxltyNMaa2xMXFERcX5+4wjLFueWOMqS179+5lzZo17g7DGEvuxhhTW1577TUSExPdHYYxltyNMcYYT+PR19yLiorIyMigsLDQ3aEY06A4HA46duxI06ZN3R2KR7HR8qah8OjknpGRQatWrQgNDbUvnTFOqkp2djYZGRl06dLF3eF4JFW1vznGrTy6W76wsJCAgAD7khlThogQEBBgPVp1wP7WmIbCo1vuYF82Yypj34u6ceONNxIdHe3uMIzx7JZ7QzBr1ix69OjB2LFj3R0K06ZNY8aMGedVx+rVqxkxYsQZ18nNzWXOnDmuzz/99BPXX3/9ee23vvzpT3/is88+q9N9VOf/4YUXXiA/P9/1efjw4eTm5tZ4X/Pnz+d3v/tdjbcz56ZHjx6MHj3aTp6M21lyr2Nz5sxh2bJlLFy4sFrrnzx5slb2q6qUlJTUSl01dXpyv/jii3n33XfdEktNPfnkkwwePNjdYVRI7suWLaN169buC8hUy549e/j3v//ttu+eMadYcq9Dd999N99//z2jRo1i5syZ5OTkcN111xEdHU1CQgJbtmwBSltykyZNYsiQIYwbN47hw4e7lsXFxfHkk08C8NhjjzF37lyOHj1KUlISvXr1Iioqin/9618ApKen06NHD37729/Sq1cvfvzxR/785z9z6aWXMnjwYHbu3FlpnEuWLCEyMpKYmBiuuOIKoHS8wsSJE4mKiiIuLo5Vq1ZV2O70FmhkZCTp6ek89NBDpKamEhsbywMPPEB6ejqRkZFnrHf+/Pn88pe/ZNiwYXTv3p2pU6ee9fjOmjWLnj17Eh0dza9+9SsAjh07xu233058fDxxcXGuYzN//nyuu+46Ro4cSZcuXfjrX//K888/T1xcHAkJCeTk5AAwYcKESk9EznbM77rrLiIiIhgyZAgFBQUAvPrqq8THxxMTE8OYMWPKJWuA1NRUevXq5fq8e/duevfuzaxZs/jpp58YNGgQgwYNAiA0NJSsrCwA3njjDaKjo4mJieG2224D4MMPP6Rfv37ExcUxePBgDh48eNbjZ2rfokWLGD58uCV3436q6hGv3r176+m+/fbbcp8HDBhQ4TV79mxVVT127Fily+fNm6eqqpmZmRWWVUfnzp01MzNTVVV/97vf6bRp01RVdeXKlRoTE6Oqqo8//rj26tVL8/PzVVX1mWee0b/+9a+al5enffr00SFDhqiq6sCBA/W7777ToqIizcvLc8UVFhamJSUlmpaWpiKia9euVVXV5ORkjYyM1GPHjmleXp6GhYXpc889VyHGyMhIzcjIUFXVn3/+WVVVZ8yYoRMmTFBV1R07dmhISIgWFBToqlWr9JprrnHFXba+iIgITUtL07S0NI2IiHCVl/1cVb3z5s3TLl26aG5urhYUFGinTp107969Zzy2QUFBWlhYWC7uhx9+WN98801XWffu3fXo0aM6b948DQsL08OHD+uhQ4fUz89PX375ZVVVnTJlis6cOVNVVcePH69LliypsK8zHfMmTZpoSkqKqqrecMMNrv1nZWW5tn/00Ud11qxZFY7bwIEDXds+/PDDrnXK/t6U/bxt2za95JJLXMuys7NVVTUnJ0dLSkpUVfXVV1/V++67T1VV582bp/fcc0+lx+/074c5f0899ZQCeuLECXeHYi4QQLJWkhOt5V6P1qxZ42ppXXnllWRnZ5OXlwfAqFGj8PHxASAxMZEvvviCNWvWcM0113D06FHy8/NJT0/n0ksvRVV55JFHiI6OZvDgwezbt8/VUuvcuTMJCQkAfPnll4wePRpfX1/8/PwYNWpUpXFdfvnlTJgwgVdffZXi4uIKsYaHh9O5c2d27dpVq8fg9HqTkpLw9/fH4XDQs2dPfvjhhzPWFR0dzdixY3nrrbfw9i4dG7pixQqmT59ObGwsAwcOpLCwkL179wIwaNAgWrVqRbt27fD392fkyJEAREVFkZ6efsZ9nemYd+nShdjYWAB69+7tqmvbtm0kJiYSFRXFwoUL2b59e4V677zzTubNm0dxcTH/+Mc/uOWWW84Yx+eff871119PYGAgAG3btgVKb/scOnQoUVFRPPfcc5Xuy9Q9u9ZuGgqPHy1f1urVq6tc5uvre8blgYGBZ1xeHVrJPM+n/hi0aNHCVRYfH09ycjJdu3blqquuIisri1dffZXevXsDsHDhQjIzM9m4cSNNmzYlNDTUdVtT2XrK1n8mf/vb31i/fj0ff/wxsbGxbN68uVpzUnt7e5frfqzOrVVnqrd58+au902aNDnr+IOPP/6YL774gg8++ICnnnqK7du3o6osXbq0wsxc69evL1e/l5eX67OXl1eFfa1fv55f//rXQOl1+JycnCqP+elxn+qWnzBhAv/85z+JiYlh/vz5lf7+jBkzhieeeIIrr7yS3r17ExAQcMafWau4f/r3v/899913H6NGjWL16tVMmzbtjPWYulWd748xdcla7vXoiiuucA2sW716NYGBgfj5+VVYr1mzZoSEhLB48WISEhJITExkxowZrmdW5+Xl0b59e5o2bcqqVauqbOFeccUVvP/++xQUFHDkyBE+/PDDStdLTU2lX79+PPnkkwQGBvLjjz+Wi3XXrl3s3bu3QsIMDQ1l06ZNAGzatIm0tDQAWrVqxZEjR856DKqq93Tjxo3j66+/LldWUlLCjz/+yKBBg3j22WfJzc3l6NGjDB06lJdeesn1xzUlJeWMdVelX79+bN68mc2bNzNq1KhqH/Oyjhw5QlBQEEVFRVUOqHQ4HAwdOpTf/OY3TJw40VVe1TFMSkpi8eLFZGdnA7jGCuTl5REcHAzAggULavzzmtphLXfTUFxQLXd3mzZtGhMnTiQ6OhpfX98z/hFOTExk5cqV+Pr6kpiYSEZGhiu5jx07lpEjR9KnTx9iY2MJDw+vtI5evXpx0003ERsbS+fOnauc0OKBBx5g9+7dqCpJSUnExMQQHh7O3XffTVRUFN7e3syfP79cCxVKW51vvPEGsbGxxMfHc8kllwAQEBDA5ZdfTmRkJFdffTX33HOPa5vf/va3Z633dFu2bCEoKKhcWXFxMbfeeit5eXmoKn/4wx9o3bo1jz32GFOmTCE6OhpVJTQ0lI8++uiM9VdHdY95WU899RT9+vWjc+fOREVFVXnCM3bsWN577z2GDBniKps0aRJXX301QUFB5QYzRkRE8OijjzJgwACaNGlCXFwc8+fPZ9q0adxwww0EBweTkJDgOtEy9evmm2+mb9++rstExriLeEr3UZ8+fTQ5Oblc2Y4dO+jRo4ebIjK14fDhw9xxxx0sWbLE3aHUmRkzZpCXl8dTTz1Vr/u170fdSE1NZd++fTgcDtfLx8eH4OBgvLy87NG0plaJyEZV7XN6uZ1emgbNz8/PoxP76NGjSU1N5fPPP3d3KKaWvPbaazzzzDMVygsKCnA4HEyePJk5c+a4kr7D4aBly5bs2LEDgP/7v/9j5cqVrmUOh4OAgACef/55AN59911SU1PLbd+mTRvXANHt27eTn59fbvsWLVrQpk2b+jsIxu0adHIXkWHAi0ATYK6qTndzSMbUqvfff9/dIZhadtdddzF48GAKCgooLCx0vZo1awbA0KFDadWqlau8oKCg3MDUoqIijhw5QmZmpquOVq1auZa/+eabfPDBB+X22bVrV1dyv/feeyucLEZHR/PNN98A8Itf/IItW7aU61VISEhwXSb89a9/zcGDB8udHMTGxroGmL7yyiucOHGi3PahoaH06VPaeNy8eTNNmzZ1be/j44Ovr+9ZL7+Z2tVgu+VFpAmwC7gKyAA2ADer6reVrW/d8sbUjH0/GqeSkpJyJw0FBQWoKt26dQNgw4YNHDx40LWssLAQf39/brzxRqD04U9paWnl6ujevbvrYVljxoxhz5495ba/8sorWbRoEQAXXXRRhYck3XLLLa5Boy1btuTYsWPllk+aNIlXXnkFVaVNmzblTgwcDgcTJkxg8uTJFBQUMHbs2HLLfHx8GDZsGIMHD+bYsWPMnz+/3ImHw+EgIiKCzp07c/z4cdLS0irU36xZM4+9FNIYu+X7AntU9XsAEVkEXAtUmtyrYte3jKmooZ7Um7Pz8vLC19cXX1/fSpfHx8efcft77733jMuXLl16xuWpqanlEn9hYSEtW7Z0LV+0aBH5+fnlTj4iIiKA0oGwEyZMqLD9qbuGjh8/zu7duyucvAQGBjJ48GAyMzMrnSvhxRdf5N5772X37t1ERUVVWD537lzuuOMONm7cyMiRI8slfofDwdNPP82gQYPYunUrzz77bIXxEuPGjaNbt26kpaWxevXqCtv36tULPz8/Dh8+TG5ubrntmzZt6pYc1JCTezDwY5nPGUC/mlTgcDjIzs62aV+NKUOd87k7HA53h2IaoRYtWlR4nkZZZ5pYytvbmxdeeKHK5a1bt2br1q1VLg8JCeHQoUMVTg5CQkIA6NixI++880655QUFBa4THn9/f0aMGFFu24KCAtfdDTk5OXz11VfllhUUFDBgwAC6devG119/ze23314hrnXr1tGvXz+WLFnCnXfeWW6Zl5cXW7dupWfPnlX+XHWhIXfL3wAMVdU7nZ9vA/qq6u/LrDMJmATQqVOn3qffe1xUVERGRobNW23MaRwOBx07dqRp06buDsWYBu1UjhQR8vPzXScXZU8gevfujZ+fH7t27WLNmjUVTj4mT57seqpkbauqW74hJ/f+wDRVHer8/DCAqlYchkrl19yNMcYYT1ZVcm/IT6jbAHQXkS4i0gz4FfDBWbYxxhhjLngN9pq7qp4Ukd8Byym9Fe51VbXZMIwxxpizaLDJHUBVlwHL3B2HMcYY05g02GvuNSUimcDZZ/OomUAgq5brvNDYMTx/dgzPnx3D82fH8PzVxTHsrKrtTi/0mOReF0QkubKBCqb67BiePzuG58+O4fmzY3j+6vMYNuQBdcYYY4w5B5bcjTHGGA9jyf3M/u7uADyAHcPzZ8fw/NkxPH92DM9fvR1Du+ZujDHGeBhruRtjjDEexpJ7JURkmIjsFJE9IvKQu+NpbEQkRERWicgOEdkuIpPdHVNjJSJNRCRFRD5ydyyNlYi0FpF3ReQ75+9kf3fH1NiIyB+c3+VtIvKOiNisQ2chIq+LyCER2VamrK2IfCoiu53/tqmr/VtyP41zHvnZwNVAT+BmEanf6Xwav5PAH1W1B5AA3GPH8JxNBna4O4hG7kXgE1UNB2Kw41kjIhIM3Av0UdVISp8Y+iv3RtUozAeGnVb2ELBSVbsDK52f64Ql94pc88ir6gng1DzypppUdb+qbnK+P0LpH9Ng90bV+IhIR+AaYK67Y2msRMQPuAJ4DUBVT6hqrluDapy8AR8R8QZ8gZ/cHE+Dp6pfADmnFV8LLHC+XwBcV1f7t+ReUWXzyFtiOkciEgrEAevdHEpj9AIwFShxcxyNWVcgE5jnvLwxV0SqnozcVKCq+4AZwF5gP5CnqivcG1Wj1UFV90NpIwhoX1c7suRekVRSZrcUnAMRaQksBaao6mF3x9OYiMgI4JCqbnR3LI2cN9ALeFlV44Bj1GFXqCdyXhe+FugCXAy0EJFb3RuVORtL7hVlACFlPnfEuqBqTESaUprYF6rqe+6OpxG6HBglIumUXhq6UkTecm9IjVIGkKGqp3qO3qU02ZvqGwykqWqmqhYB7wGXuTmmxuqgiAQBOP89VFc7suRekc0jf55ERCi9xrlDVZ93dzyNkao+rKodVTWU0t/Bz1XVWks1pKoHgB9F5FJnURLwrRtDaoz2Agki4uv8bidhgxLP1QfAeOf78cC/6mpHDXrKV3eweeRrxeXAbcBWEdnsLHvEOYWvMfXt98BC58n698BEN8fTqKjqehF5F9hE6Z0wKdjT6s5KRN4BBgKBIpIBPA5MBxaLyB2UnjTdUGf7tyfUGWOMMZ7FuuWNMcYYD2PJ3RhjjPEwltyNMcYYD2PJ3RhjjPEwltyNMcYYD2PJ3RhjjPEwltyNMcYYD2PJ3RhTbSLiIyL/cc4zH1p2ruoa1tNMRL5wzjJmjKllltyNMTVxO/CeqhafTyXO6ZRXAjfVSlTGmHIsuRtjABCRVSJylfP90yIyq5LVxlLJ87BFpKtzStV4Z4v+O+f0qttEZKGIDBaR/4rIbhHp69zsn876jDG1zJK7MeaUx4FHRWQsEAf8oexC57PZu6pq+mnll1I6A+BEVd3gLO4GvAhEA+HALcAvgPuBR5zrbAPi6+QnMeYCZ9e7jDEAqOoXzlm/7gMGVtL1HgjknlbWjtKW/JjTJlhKU9WtACKyHVipqioiW4FQ5/6KReSEiLRS1SO1/xMZc+GylrsxBgARiQKCgONVJNsCwHFaWR7wI6UzAZZ1vMz7kjKfSyjfqGgOFJ5rzMaYyllyN8YgIkHAQuBa4JiIDD19HVX9GWgiImUT/AngOmCciNxSw30GAJmqWnTOgRtjKmXJ3ZgLnIj4Au8Bf1TVHcBTwLQqVl9B6bVzF1U9BowA/iAi19Zg14OAZTUO2BhzVjafuzGm2kQkDrhPVW+rhbreAx5W1Z3nH5kxpixruRtjqk1VU4BVItLkfOpxjrz/pyV2Y+qGtdyNMcYYD2Mtd2OMMcbDWHI3xhhjPIwld2OMMcbDWHI3xhhjPIwld2OMMcbDWHI3xhhjPIwld2OMMcbD/D/Qqsths9V4swAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from os import getcwd, chdir\n",
    "cwd = getcwd()\n",
    "chdir(r'../../timml')\n",
    "import timml\n",
    "chdir(cwd)\n",
    "\n",
    "M0 = timml.Model3D(kaq=K1, z=[0.,-b])\n",
    "timml.inhomogeneity1d.StripInhom3D(M0, x1=2*L/3., x2=L, kaq=K2, z=[0,-b], topboundary='conf')\n",
    "timml.ImpLineDoublet1D(M0, xld=0.)\n",
    "timml.HeadLineSink1D(M0, xls=L, hls=0.)\n",
    "timml.StripAreaSink(M0, 0., L, R)\n",
    "M0.solve(silent=True)\n",
    "H0 = M0.headalongline(X, 0.).flatten()\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "plt.subplots(figsize=[8,2])\n",
    "plt.plot(X/1000., H0, 'k--', mfc='none', label='forward solution, semi-analytical')\n",
    "plt.xlabel('$x$ (km)')\n",
    "plt.ylabel('$h$ (m)')\n",
    "plt.legend(loc=3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgAAAACdCAYAAADCH8BOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnuUlEQVR4nO3deVxVdRr48c8DiIopKmoqJLhUKAiouJCBmuEylalpWuaI/ix1bJ3JxnIaxykby6bFFrX8pVNjqdmi/VKnbHTEyVRcc6s0VFAHEEVZElme3x8XbyCouMBled6vFy/uOd+zPOdwL+e53/M936+oKsYYY4ypXtxcHYAxxhhjyp8lAMYYY0w1ZAmAMcYYUw1ZAmCMMcZUQ5YAGGOMMdWQJQDGGGNMNeTh6gDKU6NGjTQgIMDVYRhjjDHlYsuWLcdVtXFJZZU6ARCRfsDrgDswT1VnXGz5gIAA4uLiyiU2Y4wxxtVE5NCFyirtLQARcQfeAvoD7YD7RKSda6MyxhhjKodKmwAAXYD9qvqzqp4FFgF3uzgmY4wxplKozLcAfIGEQtOJQFcXxWKMMZeUl5fH2LFjOXLkCCtXrsTd3Z033niDL774oshyNWrU4MsvvwTgpZdeYvXq1UXKvb29+fjjjwH4y1/+wrffflukvGnTprz//vsA/PGPf2Tbtm1Fylu1asWcOXMAePTRR9m3b1+R8qCgIF599VUAHnzwQQ4dKlqLHB4ezgsvvADAiBEjSElJKVIeGRnJs88+C8DgwYPJyMgoUt63b1/+8Ic/ANCvXz/y8/OLlA8cOJDf/e53nD17ljvvvJPz3X///cTExJCWlsa9995brHzMmDEMHz6cY8eOMWrUqGLlDz/8MAMGDODAgQNMmDChWPmkSZOIjo5m165d/P73vy9WPnXqVLp3787mzZuZMmVKsfIZM2bQsWNH1q1bx/PPP1+s/PXXX6dt27asWrWKV155xTm/devWzJ49u9jyZaUyJwBSwrxiAxuIyEPAQwAtWrQotkJOTg6JiYmcOXPmmgdojLk2atWqhZ+fHzVq1HB1KFclOTmZBQsW0LJlS86Nw5KdnV3sAunp6el8febMmWLlHh4eFy3Pyspyvv7ll18uu/yXX34psuzllhf+f5qZmXnJ8ry8vCLl2dnZAKhqsXVLU56TkwNAfn7+FZXn5uYCjoTtSsrPHU9ubu5llRf+u5QHqayDAYlIBPAXVe1bMP00gKr+7ULrhIeH6/mNAOPj46lbty4+Pj6IlJRTGGNcSVVJTU0lPT2dli1bujqcq3Ls2DGaN2/OnDlzGDdunKvDMdWAiGxR1fCSyipzG4DNwI0i0lJEPIHhwPLL3ciZM2fs4m9MBSYi+Pj4VIlausr6hctUTZX2FoCq5orIw8C/cDwG+J6q7r6SbdnF35iKrap9Rqva8ZjKqTLXAKCqK1T1JlVtrarTXR2Pq3z++efs2bPnstdbu3ZtiQ1sLua1114rcp/qN7/5DWlpaZe97/I0Z84cZ4OoBQsWcPToUWfZ2LFjL3nubrnlFgAOHjzIhx9+WHaBXoGjR48yZMiQMt/Pddddd9Hy7du3s2LFCuf08uXLmTHjot1yXFBAQADHjx+/onUrukaNGrFx40YGDhzo6lCMqdwJgHE0IrmSBOBcI5bLdX4CsGLFCurXr39F2yov48eP57e//S1QPAGYN28e7dpdvPuIcy2sK2IC0Lx5c5YuXerqMIolAAMGDGDy5MkujKhi8vT0pEuXLjRp0sTVoRhjCYCrHTx4kMDAQEaNGkVISAhDhgxxXmC3bNlCjx496NSpE3379uXYsWMA9OzZk2eeeYYePXrw4osvsnz5ciZNmkRYWBgHDhygZ8+ezh4Pjx8/zrnujxcsWMDQoUO566676NOnDwCnT59m0KBBtGvXjvHjxzsfx5kwYQLh4eEEBQUxdepUAGbNmsXRo0fp1asXvXr1Aop+W3vllVcIDg4mODiY1157zXl8bdu25cEHHyQoKIg+ffoUaUF8zscff0xwcDChoaFERUUBjpaykyZNonPnzoSEhDB37lzAUXPRs2dPhgwZQmBgICNGjHDeW508eTLt2rUjJCSEJ598EnA8JvXyyy+zdOlS4uLiGDFiBGFhYfzyyy/OczV79myeeuopZzwLFizgkUceAX799jt58mRiY2MJCwvj1VdfJTIyku3btzvX6d69Ozt37rzg3/rYsWNERUURFhZGcHAwsbGxAHz11VdERETQsWNHhg4d6mwVHBAQwDPPPENERATh4eFs3bqVvn370rp1a+cjXAcPHiQ4OLjE/b377rt07tyZ0NBQ7rnnHuf7KiYmhkcffZRbbrmFVq1aOROIjIwMevfuTceOHWnfvj3Lli0rts2RI0cWmT9ixAiWL1/On//8ZxYvXkxYWBiLFy9mwYIFPPzwwwAkJSUxaNAgQkNDCQ0NdSZUAwcOpFOnTgQFBfHOO+9c8LxVJZmZmcybN6/YY3fGuISqVpufTp066fn27NlTbF55io+PV0DXr1+vqqqjR4/WmTNn6tmzZzUiIkKTk5NVVXXRokU6evRoVVXt0aOHTpgwwbmNUaNG6ccff+yc7tGjh27evFlVVVNSUtTf319VVefPn6++vr6ampqqqqpr1qzRmjVr6oEDBzQ3N1dvv/1253bOLZObm6s9evTQHTt2qKqqv7+/pqSkOPd1bjouLk6Dg4M1IyND09PTtV27drp161aNj49Xd3d33bZtm6qqDh06VD/44INi5yE4OFgTExNVVfXkyZOqqjp37lx97rnnVFX1zJkz2qlTJ/355591zZo1Wq9ePU1ISNC8vDzt1q2bxsbGampqqt50002an59fZDtTp07VmTNnFjs3haeTk5O1devWzvn9+vXT2NhYVVWtU6eO83zdcccdzmUWLFigjz32mKqq/vDDD1rS+6uwl19+WZ9//nnneT19+rSmpKRoZGSkZmRkqKrqjBkzdNq0ac5z+/bbb6uq6uOPP67t27fX06dPa3JysjZu3FhVHe+foKCgEvd3/Phx5+spU6borFmzVNXxfhkyZIjm5eXp7t27ncedk5Ojp06dUlXH+6Z169bOc3nuHKxdu1bvvvtuVVVNS0vTgIAAzcnJ0fnz5+vEiROd+ys8fe+99+qrr77qPO60tDRV/fU9lpWVpUFBQc54z3+PnePqz+q1kJCQoIC+++67rg7FVBNAnF7gmlhpGwGWlZ49exabd++99/K73/2OrKwsfvOb3xQrj4mJISYmhuPHjxe7H7t27dpL7vOGG26ge/fuADzwwAPMmjWLfv36sWvXLqKjowHHt+FmzZo51xk2bNhlHNWvoqOjadiwoXO6S5cutGrVCoD77ruP9evXM2TIEJYsWcI777xDbm4ux44dY8+ePYSEhFxwu+vXr2fQoEHUqVMHcHT+ERsby4ABA2jZsiVhYWEAdOrUiYMHDxZbv3v37sTExHDvvfcyePBgwPHNeOfOnc5vqKdOneKnn35yVqP6+fkBEBYWxsGDB+nWrRu1atVi7Nix3HHHHZfVvqFx48a0atWK7777jhtvvJEffvjB+Te5kKFDh/Lcc88xc+ZM3nvvPWJiYi66fOfOnRkzZgw5OTkMHDiQsLAw/vOf/7Bnzx7nvs6ePUtERIRznQEDBgDQvn17MjIyqFu3LnXr1qVWrVqXbHuxa9cu/vSnP5GWlkZGRgZ9+/Z1lg0cOBA3NzfatWtHUlIS4Pgy8Mwzz7Bu3Trc3Nw4cuQISUlJNG3a1Llejx49mDhxIsnJyXz66afcc889RZ5JL8m///1vZxsMd3d3vL29AUeN0meffQZAQkICP/30Ez4+PhfdVmWn9hSAqUAsAagAzm8RLCKoKkFBQWzYsKHEdc5daEvi4eHhrMo//9Gp89crad/x8fG8/PLLbN68mQYNGhATE3PJR7Au9o+tZs2aztfu7u4l3gKYM2cOGzdu5MsvvyQsLIzt27ejqrzxxhtFLlzgSKrO32Zubi4eHh5s2rSJb775hkWLFvHmm2/y73//+6JxFzZs2DCWLFlCYGAggwYNumRLbS8vL6Kjo1m2bBlLliy55EBTUVFRrFu3ji+//JKRI0cyadIkGjRoQHR0NB999FGJ65w7Tjc3tyLH7ObmVqwdx+jRo9m2bRvNmzdnxYoVxMTE8PnnnxMaGsqCBQuKJKOFt3Xub7dw4UJSUlLYsmULNWrUICAgoMS/+8iRI1m4cCGLFi3ivffeu+gxX8jatWtZvXo1GzZswMvLi549e1aJx/xKy54CMBWBJQDnudg3di8vr4uWN2rUqFTf+M93+PBhNmzYQEREBB999BG33norN998MykpKc75OTk5/PjjjwQFBRVbv27duqSnpzunAwIC2LJlC126dLlkA7FNmzYRHx+Pv78/ixcv5qGHHuL06dPUqVMHb29vkpKSWLlypbNm5Ny+GjVqVGQ7UVFRxMTEMHnyZFSVzz77jA8++KDU5+DAgQN07dqVrl278sUXX5CQkEDfvn2ZPXs2t912GzVq1ODHH3/E19f3gtvIyMhw1tJ069aNNm3aXPJcFTZ48GCmT5+Ov78/L774YqnWHTt2LHfddReRkZHOmpVNmzbx5ptvOr/1nnPo0CF8fX158MEHyczMZOvWrUyZMoWJEyeyf/9+2rRpQ1ZWFomJidx0002XPGfnmz9/fpHp9PR0mjVrRk5ODgsXLrzouQNHDUuTJk2oUaMGa9asKdb96zkxMTF06dKFpk2bOt+PFzuvvXv3Zvbs2Tz++OPk5eWRmZnJqVOnaNCgAV5eXuzbt4/vvvvuso+3MrIaAFORWCPACqBt27b84x//ICQkhBMnTjBhwgQ8PT1ZunQpf/zjHwkNDSUsLKxYf9/nDB8+nJkzZ9KhQwcOHDjAk08+yezZs7nlllsu+ThVREQEkydPJjg4mJYtWzoba3Xo0IGgoCDGjBlTpCr8oYceon///s5GgOd07NjReWHo2rUrY8eOpUOHDqU+B5MmTaJ9+/YEBwcTFRVFaGgoY8eOpV27dnTs2JHg4GDGjRt30acX0tPTufPOOwkJCaFHjx7OvswLi4mJYfz48c5GgIU1aNCAdu3acejQIbp06VJs3ZCQEDw8PAgNDXVuu1OnTtSrV4/Ro0c7lzt8+DC1a9cutv7atWsJCwujQ4cOfPLJJzz22GM0btyYBQsWcN999xESEkK3bt2uWQOx5557jq5duxIdHU1gYOAllx8xYgRxcXGEh4ezcOHCC65z/fXX07Zt2yLH3KtXL/bs2eNsBFjY66+/zpo1a2jfvj2dOnVi9+7d9OvXj9zcXEJCQnj22Wfp1q3b1R1sJWM1AKYiqLRdAV+JkroC3rt3L23btnVRRI5W3HfeeSe7du1yWQzmyh09epSePXuyb98+3Nwc+fSkSZMYOXLkRdtMVGZZWVm0b9+erVu3Ou/nlwdXf1avhZycHPbv30+zZs0q/OOzpmqoql0BG+NS77//Pl27dmX69OnOiz/AzJkzq+zFf/Xq1QQGBvLII4+U68W/qqhRowZt27a1i7+pEKwGoAp8qzCmOqgKn9VTp04xf/58+vTpc8kOqIy5FqwGwBhjKoATJ07wxBNPsHnzZleHYowlAGAtc42p6KrKZ7SqHIepGqp9AlCrVi1SU1Ptg2lMBaWqpKamUqtWLVeHcs3YUwCmIqj2/QD4+fmRmJhISkqKq0MxxlxArVq1nD0/Vmb2RcNUJNU+AahRowYtW7Z0dRjGmGrEagBMRVDtEwBjjCkv/v7+JCQk0KBBA1eHYowlAMYYU148PDyqxK0MUzVU+0aAxhhTXk6cOMFzzz3H999/7+pQjLEEwBhjyktqaip//vOf2blzp6tDMcYSAGOMKS/2FICpSCwBMMaYcmZPAZiKwBIAY4wpJ1YDYCoSSwCMMaacWQ2AqQjsMUBjjCknbdq04eTJk3h5ebk6FGMsATDGmPLi7u5O/fr1XR2GMYDdAjDGmHKTlJTEU089xY4dO1wdijGWABhjTHlJTU1l5syZ/PDDD64OxRhLAIwxprzYUwCmIrEEwBhjypk9BWAqAksAjDGmnFgNgKlILAEwxphyZjUApiK4ZAIgIu5lGYCIzBSRfSKyU0Q+E5H6hcqeFpH9IvKDiPQtNL+TiHxfUDZL7NNkjKkEgoKCyM3NZfDgwa4OxZhS1QDMK+MYvgaCVTUE+BF4GkBE2gHDgSCgH/B2oWRkNvAQcGPBT78yjtEYY66aiODu7o6bm1W+GtcrzbvQXUT+ev5MEakrIh9fbQCq+pWq5hZMfgf4Fby+G1ikqtmqGg/sB7qISDOgnqpuUMcNtfeBgVcbhzHGlLUjR44wYcIEtm3b5upQjClVAjAGiBCRsedmiEgoEAfsvcbxjAFWFrz2BRIKlSUWzPMteH3+fGOMqdBSU1OZM2cO8fHxrg7FmEt3BayquSIyGFgrIkdxXGynAGNVdXVpdiIiq4GmJRRNUdVlBctMAXKBhedWKymci8y/0L4fwnG7gBYtWpQmXGOMKRP2FICpSC6ZAIjI28BO4CVgPvA90FVVk0q7E1W9/RL7GAXcCfTWXz8hicANhRbzA44WzPcrYf6F9v0O8A5AeHi4ffqMMS5n7ZZNRVCaWwDbgWDgd0AN4CbgXRGZLiLDrzYAEekH/BEYoKpZhYqWA8NFpKaItMTR2G+Tqh4D0kWkW0Hr/98Cy642DmOMKWtWA2AqktLcAnin8LSI+AEhQHvgDmDRVcbwJlAT+LogK/5OVcer6m4RWQLswXFrYKKq5hWsMwFYANTG0WZgZbGtGmNMBSMi1K5dG3f3Mn262phSkeqUkYaHh2tcXJyrwzDGGGPKhYhsUdXwksouWQNQsIGawFNABJAJbAWWqOqBaxalMcYYY8pNaXujmAn44+ikJxHwAeaLyANlFZgxxlQ18fHxPPDAA2zdutXVoRhT6gQgFHhSVXcAXVT1SSAamFhmkRljTBVz4sQJFi5cyJEjR1wdijGluwWAowbgXLJwXEReAfYB+WUSlTHGVEHVqc2VqfhKMxhQBPClqp4omDUURy+A1xW8NsYYcxmsHwBTEZSmBmAU8JaI/AisAlap6odlG5YxxlQ9VgNgKpLS9AMwHkBEAoH+wAIR8QbW4EgI/lvo+XxjjDEX4OHhQZMmTahZs6arQzHmyvoBEJHaQC8cCUHEhZ4xrGisHwBjjDHVyVX3A1CwkduAEUAasAvH+ABPqmr2tQjSGGOMMeWntI8BAvwT+H/Ad0Ar4M/A7rIIyhhjqhJVJS4ujtGjR9O7d2+2bNni6pCMKX0NALBfVT8reP1xWQRjjDFVyaFDh3jllVf4/PPPOXz4MO7u7vTq1Ysbbrjh0isbU8ZK8xjg+yLyOLBBRP5Q9iEZY0zllJ2dzYoVK9i4cSMAubm5zJ07l9DQUObPn09SUhJff/01TZo0cXGkxpSuBuAfOHoCbAr0FZFHgR0FPztV1WoDjDHVVlpaGqtWrWLZsmV8+eWXpKenM2LECLp27Urr1q05ceIEXl5erg7TmGJK8xjgN8A356ZFxANohyMp6IrdDjDGVDOpqan4+PgAEBUVxffff0/jxo0ZNmwYgwcP5rbbbnMuaxd/U1Fd7WiAH5RhbMYYUyHk5+ezefNmli9fzvLly0lISCA5ORlPT09eeuklvL296dKlC+7u7q4O1ZhSu5yxALxwjAYYw6+jAc5V1YVlFJsxxrjc0qVLefjhh0lKSsLd3Z3IyEjGjBlDTk4Onp6e9OvXz9UhGnNFSpsAhAJ3q2qaiHRR1e4FtQJrAUsAjDFVwtGjR1mxYgXLly/nySefJCoqCj8/P3r06MGAAQPo378/DRs2dHWYxlwTNhqgMaZay8jI4IUXXmDFihXs2LEDAH9/f1JSUgDo1q0bixcvdmWIxpSJy+4KWEQ8gSE4ngpYpKpHyyKwsmBdARtj/ve//7Fq1SpEhFGjRpGXl4evry+BgYH079+f/v370759exuxz1QJ16Qr4HNU9SxgowEaYyqNuLg4li1bxooVK9i6dSsAt956K6NGjcLd3Z1Dhw7ZAD2m2rmcroCNMaZS+N///seSJUucw+/OmjWLF154gdq1azN9+nS2bdvGunXrnMvbxd9UR1c0GmBlZbcAjKmasrKyiI2N5euvv+brr79m586dAOzdu5fAwEAOHz7MddddZw34TLVzTW8BGGOMq+Xn57N9+3aaNm1K8+bNWbFiBUOHDsXT05PIyEhmzJhBdHQ0N910EwAtWrRwccTGVDyWABhjKoWEhATnN/zVq1dz/PhxXnrpJSZNmkR0dDSrVq0iMjLSet4zppQsATDGVEgnT54kKSmJwMBAsrKyaN26NTk5OTRt2pT+/fsTHR1Nnz59APD29qZv374ujtiYysUSAGNMhXDq1CliY2NZs2YNa9asYfv27XTv3p3Y2Fi8vLx4//33CQoKIjg42B7RM+YasKcAjDFlZ+FCCAgANzfH74W/dhyakZFBbGysc/r+++/nrrvu4q233qJ+/fpMmzaNF1980Vk+fPhwez7fmGvIngIwxpSNhQvhoYcgK8s5K69mTT7p25fXUlLYvHkzubm5pKSk0KhRI9avX09OTg4RERHUqlXLhYEbU3XYUwDGmHKXN3ky7oUu/gDu2dl0Xb4c7daNSZMm0atXL+rWrQs4OuYxxpQfSwCMMddEeno6X3zxBbGxscTGxrIzMbHE5VqIsGHDhnKOzhhzvgrTBkBEnhQRFZFGheY9LSL7ReQHEelbaH4nEfm+oGyW2E1BY8qVqrJnzx7mzp3Lv/71L8BxT3/EiBEsXLgQPz8/Tnt7l7iu2DP5xlQIFSIBEJEbgGjgcKF57YDhQBDQD3hbRNwLimcDDwE3FvzYgNzGlIOXX36ZAQMG0LhxY4KCghg/fjwfffQRAM2aNWPnzp2cPHmSVatWUf+tt+D8Z/K9vGD6dBdEbow5X0W5BfAq8BSwrNC8u3GMNpgNxIvIfqCLiBwE6qnqBgAReR8YCKws14iNqaJUlZ9++onvvvuODRs2kJ+fz9y5cwFYvHgx6enp3HXXXURGRhIZGUmbNm2c67Zv3/7XDY0Y4fg9ZQocPgwtWjgu/ufmG2NcyuUJgIgMAI6o6o7zavJ9ge8KTScWzMspeH3+fGPMFcjMzKROnToATJ06lTfffJMTJ04AUK9ePW677TbnsuvXr7+8gXNGjLALvjEVVLkkACKyGmhaQtEU4BmgT0mrlTBPLzL/Qvt+CMftAusP3FR7Z8+eZefOncTFxREXF8emTZvYu3cvx48fx9vbm+bNmzNw4EAiIiKIiIigbdu2uLn9eqfQRs0zpuoolwRAVW8vab6ItAdaAue+/fsBW0WkC45v9jcUWtwPOFow36+E+Rfa9zvAO+DoB+DKj8KYyiU3N5c9e/YQFxdH//79adasGfPmzWPixIkA+Pj40LlzZ+655x5yc3MBGDduHOPGjXNl2MaYcuLSWwCq+j3Q5Nx0wf39cFU9LiLLgQ9F5BWgOY7GfptUNU9E0kWkG7AR+C3wRvlHb0zFc/jwYf7+978TFxfHtm3b+OWXXwBYtGgRw4YN44477qBJkyaEh4fj7+9vveoZU425vA3AhajqbhFZAuwBcoGJqppXUDwBWADUxtH4zxoAmmojOzubXbt2sX37dufPAw88wLhx41BV5s2bR8eOHRk3bhydO3cmPDzc2VDP398ff39/Fx+BMaYiqFAJgKoGnDc9HSj2zJCqxgHB5RSWMS6TmprKjh078PDwICoqiuzsbLy9vcnOzgbguuuuIywszNmIr0WLFpw6dQoPjwr10TbGVED2X8KYCkBVndXxf//731m3bh3btm0jISEBgD59+hAVFUXNmjX529/+xg033EBYWBitWrUq0khPROzib4wpFftPYUw5O3LkCNu3b2f37t3s2rWLXbt24e7uzubNmwH46quvSEhI4NZbbyUsLIwOHToQFhbmXP+JJ55wUeTGmKrEEgBjyoCqkpyc7LzA//TTT7zxxhuICE8//TQffPABAH5+fgQHBxe5wK9cubLIt3pjjCkLNhywMVchJyeH+Ph49u3bx+23346Xlxdz5szhmWee4eTJk87lfHx82LdvH40aNWL79u1kZmYSFBRE/fr1XRe8MabKs+GAjblKaWlpeHp64uXlxcaNG5kxYwb79u1j//79zmfoN23aROfOnWndujXDhg0jMDCQ4OBggoODadKkifMef+Fv+8YY4yqWABiDo8o+Pz8fd3d3jh07xrvvvsuBAwfYv38/+/fvJzk5mU8//ZRBgwaRnZ3Njz/+SNu2bRk0aBCBgYHOiz1AdHQ00dHRLj4iY4y5OEsATLWTlZXF4sWL2b9/v/Mif+DAAaZNm8ajjz5KRkYGU6dOxc/PjzZt2nDXXXcVucBHRUWxe/duFx+FMcZcHUsATJWjqnz44YccOnSIgwcPOn8PGTKE6QVD0Y4ZMwZ3d3cCAgJo06YNXbt2dY5k17p1a7Kysqhdu7YrD8MYY8qUJQCm0sjOznYORrNkyRL27t3rvMAfOnSIiIgI/vnPfyIiPPbYY6SmptKkSRP8/f0JCQnh5ptvBsDLy4uff/4ZPz8/atSoUWw/bm5udvE3xlR5lgAYl8vLyyM5OZm0tDTatm0LwNtvv83GjRs5cuQIiYmJHDlyhKCgIL77zjFC9N/+9jd27NhBs2bN8Pf3p2vXrnTv3t25zbi4OJo0aYKXl1eJ+2zZsmXZH5gxxlRglgCYMpGfn8+JEydISkpy/pw+fdo50tz06dP57LPPSEpK4tixY+Tl5eHr60tiYiIAq1evJi4uDl9fX4KDg+nXrx/t2rVzbn/lypU0aNDggsPTBgQElPkxGmNMZWYJgLkkVSUzM5PU1FR8fX3x8PBgx44d/Pe//yU1NZXU1FROnDhBcnIyy5cvx9PTk8cff5w33ig6SKObmxtjx47F3d0dNzc3GjduTHBwMM2bN8fPz48WLVo4l/3kk08uOlJd06ZNy+x4jTGmOrAEoBpQVTIyMjh9+jQNGzakdu3aJCYm8u2333L69GlOnTrFqVOnOH36NJMmTcLX15ePP/6YadOmOS/uZ8+eBeDnn3+mZcuWrFy5kqeffhqAunXr4uPjQ+PGjUlPT8fHx4dBgwbRpk0brr/+epo0acL111/P9ddf7+zh7umnn3auXxIbptYYY8qWJQAucq6K/OzZs5w9e5bs7GzOnj1L06ZNady4MRkZGXz77bfO8jNnzpCZmUlkZCSBgYHEx8cza9YssrKyyMrKIjMzk6ysLP70pz9x66238p///IeRI0dy+vRp0tPTyc/PBxxV67179+bbb79l2LBhznhEhLp16zJy5Eh8fX3x9vbmpptuwsfHBx8fHxo2bIiPjw8NGjQAYPz48cTExNCwYUM8PT2LHV+vXr3o1atX+ZxMY4wxl80SgCuUm5tLixYtyM/PL/LzyCOPMG3aNNLS0mjWrFmx8r/+9a88++yzJCYmljgu+2uvvcZjjz3G4cOH6du3b7Hyd999l8DAQFJTU5k3bx516tTBy8vL+fvcMLGNGjWid+/e1K1bF29vb7y9valXr56zJXyfPn3YtWsX9erVw9vbm+uuu65I//N9+vShT58+Fzx+68LWGGMqN0sArpCbmxt33nknbm5uRX46duwIQM2aNXn00UeLlUdGRgLQsGFDZs2ahaenJzVr1sTT0xNPT09nN7EBAQGsX7/eWVazZk3q1KmDj48PAOHh4aSnp18wvqCgIObPn3/B8vr169tF3BhjqjEbDMgYY4ypoi42GJCNOWqMMcZUQ5YAGGOMMdWQJQDGGGNMNVSt2gCISApw6BpushFw/Bpurzqyc3j17BxeG3Yer56dw6t3rc+hv6o2LqmgWiUA15qIxF2ocYUpHTuHV8/O4bVh5/Hq2Tm8euV5Du0WgDHGGFMNWQJgjDHGVEOWAFydd1wdQBVg5/Dq2Tm8Nuw8Xj07h1ev3M6htQEwxhhjqiGrATDGGGOqIUsAroCI9BORH0Rkv4hMdnU8lZGI3CAia0Rkr4jsFpHHXB1TZSUi7iKyTUT+n6tjqYxEpL6ILBWRfQXvxwhXx1TZiMgTBZ/jXSLykYjUcnVMlYGIvCciySKyq9C8hiLytYj8VPC7QVnt3xKAyyQi7sBbQH+gHXCfiLRzbVSVUi7wB1VtC3QDJtp5vGKPAXtdHUQl9jqwSlUDgVDsXF4WEfEFHgXCVTUYcAeGuzaqSmMB0O+8eZOBb1T1RuCbgukyYQnA5esC7FfVn1X1LLAIuNvFMVU6qnpMVbcWvE7H8U/X17VRVT4i4gfcAcxzdSyVkYjUA6KA/wugqmdVNc2lQVVOHkBtEfEAvICjLo6nUlDVdcCJ82bfDfyj4PU/gIFltX9LAC6fL5BQaDoRu3BdFREJADoAG10cSmX0GvAUkO/iOCqrVkAKML/gNso8Eanj6qAqE1U9ArwMHAaOAadU9SvXRlWpXa+qx8DxRQloUlY7sgTg8kkJ8+xRiiskItcBnwCPq+ppV8dTmYjInUCyqm5xdSyVmAfQEZitqh2ATMqwyrUqKrhHfTfQEmgO1BGRB1wblSkNSwAuXyJwQ6FpP6y664qISA0cF/+Fqvqpq+OphLoDA0TkII5bUbeJyD9dG1Klkwgkquq52qelOBICU3q3A/GqmqKqOcCnwC0ujqkySxKRZgAFv5PLakeWAFy+zcCNItJSRDxxNHZZ7uKYKh0RERz3Xfeq6iuujqcyUtWnVdVPVQNwvA//rar2zesyqOr/gAQRublgVm9gjwtDqowOA91ExKvgc90ba0h5NZYDowpejwKWldWOPMpqw1WVquaKyMPAv3C0dn1PVXe7OKzKqDswEvheRLYXzHtGVVe4LiRTTT0CLCxI6H8GRrs4nkpFVTeKyFJgK46ne7ZhPQKWioh8BPQEGolIIjAVmAEsEZH/gyO5Glpm+7eeAI0xxpjqx24BGGOMMdWQJQDGGGNMNWQJgDHGGFMNWQJgjDHGVEOWABhjjDHVkCUAxhhjTDVkCYAxxhhTDVkCYIy5ZkSktoj8R0TcRSSg8Djnl7kdTxFZVzC6nDGmDFgCYIy5lsYAn6pq3tVspGCo7W+AYdckKmNMMZYAGGNKRUTWiEh0wevnRWRWCYuNoIS+y0WkVcFwu50Lagb2FQy9u0tEForI7SLyXxH5SUS6FKz2ecH2jDFlwBIAY0xpTQWmiMgIoAPwROHCgr70W6nqwfPm34xj1MfRqrq5YHYb4HUgBAgE7gduBZ4EnilYZhfQuUyOxBhjgwEZY0pHVdcVjPb2e6BnCdX8jYC08+Y1xlEjcM95g2bFq+r3ACKyG/hGVVVEvgcCCvaXJyJnRaSuqqZf+yMypnqzGgBjTKmISHugGZB9gQvyL0Ct8+adAhJwjP5YWHah1/mFpvMp+sWkJnDmSmM2xlyYJQDGmEsSkWbAQuBuIFNE+p6/jKqeBNxFpHAScBYYCPxWRO6/zH36ACmqmnPFgRtjLsgSAGPMRYmIF/Ap8AdV3Qs8B/zlAot/heNevpOqZgJ3Ak+IyN2XsetewIrLDtgYUyqiqq6OwRhTRYhIB+D3qjryGmzrU+BpVf3h6iMzxpzPagCMMdeMqm4D1oiI+9Vsp+CJgs/t4m9M2bEaAGOMMaYashoAY4wxphqyBMAYY4yphiwBMMYYY6ohSwCMMcaYasgSAGOMMaYasgTAGGOMqYYsATDGGGOqof8PO+so6ExWOsUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "M1 = timml.Model3D(kaq=K1+K1*dpar, z=[0., -b])\n",
    "timml.inhomogeneity1d.StripInhom3D(M1, x1=2*L/3., x2=L, kaq=K2, z=[0,-b], topboundary='conf')\n",
    "timml.ImpLineDoublet1D(M1, xld=0.)\n",
    "timml.HeadLineSink1D(M1, xls=L, hls=0.)\n",
    "timml.StripAreaSink(M1, 0., L, R)\n",
    "M1.solve(silent=True)\n",
    "H1 = M1.headalongline(X, 0.).flatten()\n",
    "dhdK = (H1-H0)/(K1*dpar)\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "plt.subplots(figsize=[8,2])\n",
    "plt.plot(X/1000., dhdK, 'k--', mfc='none', label='perturbation sensitivity, semi-analytical')\n",
    "plt.plot(ocol/1000., dhdK[ocol], 'ro')\n",
    "plt.xlabel('$x$ (km)')\n",
    "plt.ylabel(r'$\\partial h / \\partial K$')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src onerror=\"var cell = this.closest('.code_cell');var editor = cell.querySelector('.input_area');editor.style.background='rgba(0, 200, 0, 0.2)';this.parentNode.removeChild(this)\">"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-374.962504\n"
     ]
    }
   ],
   "source": [
    "set_background('rgba(0, 200, 0, 0.2)')\n",
    "print('%.6f'% dhdK[ocol])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2c. Numerical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAACdCAYAAACpbTGvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmMUlEQVR4nO3deVyVZfr48c/FLhmggIpogomgP1BEFBM1yzUqNbRs3Je2qWy+Wab100lbJsevM9oytkyUlk622LiUuaapuSWu5L6gkguLS4qmgPfvDw7nh4ILsjyc4/V+vc6Lc+5nu87Tseu5n+dexBiDUkoppZyHi9UBKKWUUqpsaXJXSimlnIwmd6WUUsrJaHJXSimlnIwmd6WUUsrJaHJXSimlnIyb1QGUlYCAABMSEmJ1GEoppVSFSU5OzjTGBF5Z7jTJPSQkhA0bNlgdhlJKKVVhRORgceV6W14ppZRyMprclVJKKSfjNLfly9rChQuZNGkSLi4uuLq62l+TJ0+mbt26LFq0iBkzZuDm5kaVKlXsr+eff55q1aqxZcsWtmzZYi+vWrUqvr6+REVF4ebmhjEGEbH6ayqllHJCmtyv4sKFC5w8eZJLly6Rl5dnf128eBGA3377jZ9++omcnBz++OMPzp8/z/nz53nyySepVq0as2fPZuzYsUX2e+rUKXx9fRk1ahTvvvsuvr6++Pr64uPjg5+fH/Pnz8fNzY3Fixezb98+AgMDCQgIIDAw0P5SSimlrkWcZeKY2NhYY3WDuoJzKSKcPn2azMxMe9I/e/Ysp0+fplu3bri4uDB//nyWLVvG6dOn7a8//viD5cuXA9C/f3+mT59+2f79/Pw4efIkAK+88gpbtmwhODiY2rVrExwcTGhoKB07dqzQ76yUUso6IpJsjIm9slxr7mWo8G32ghr51SQkJJCQkHDV5UlJSUyYMIGMjAz7q+CuAeRfSBw9epTk5GTS09MxxhAZGcm2bdsA6NChA6mpqYSGhhIaGkpISAjNmjWzH1MfCyillPPSmrsTyMnJ4ejRo2RnZ9OoUSMA/vd//5fk5GRSU1M5cOAA6enp3HfffcyfPx+Axo0bY4whNDSUsLAwwsPDiYuLo3nz5lZ+FaWUUiWgNXcn5u7uzh133HFZ2YgRIy77fO7cOc6cOWP/3K1bN/bu3cv+/ftZsWIF2dnZPPHEE3z44Yfk5eVx1113ERoaSsOGDQkPDyciIoLGjRvj7e1dId9JKaXUzdPkfovw9va+LDGPHz/e/t4Yw5EjR7h06RIAZ86cITAwkI0bNzJr1izy8vIAGDduHH/96185efIkkyZNIjIykqioKMLCwnBz05+SUkpVFvp/ZIWIEBwcbP/s5+fH999/D8DFixfZv38/27dvt9/y3717N2+++ab9YsDDw4NGjRoxefJk2rdvT3Z2NhcvXqRatWoV/2WUUkppclfX5uHhQUREBBEREfayuLg4srOz2blzJ9u2bSMlJYVt27bh5+cHwOzZs+nXrx/169cnJiaG5s2bExMTQ9u2balSpYpF30QppW4d2qBOlbmdO3cye/ZskpOT2bhxI/v37wfg8OHD1KlTh7lz57J9+3buuusuYmNjue222yyOWCmlHJM2qFMVJiIiglGjRtk/nzx50t4nH2DZsmVMnjwZAFdXV5o0aUJ8fDzvvPOOds9TSqkyoDV3ZYmsrCzWrVvHmjVrWLt2LRcvXuSnn34CoHfv3mRnZ3PXXXfRrl07WrRogZeXl8URK6VU5aM1d1Wp+Pv7X3Ugn6CgIBYvXmxv1Ofp6cmf//xnJk2aBOQPDezp6Vmh8SqllCPR5K4qnYJb9idOnGDlypWsWLHC3lL/zJkz1KxZk+joaNq1a0e7du1o06YNPj4+FkaslFKVi96WVw4lKyuLiRMnsmLFCn755RdycnJwdXVl2rRp9O3b1z5Er4eHh8WRKqVU+dPb8sop+Pv789ZbbwH5o+6tXbuWH3/8kdjY/N/27NmzGTp0KHfffTedOnWiY8eONG7cWBvqKaVuKVpzV05l48aNJCUlsXjxYvbs2QPkP8PfsmULgYGBXLp0CRcXF4ujVEqpsmFpzV1E6gKfAbWAS8BHxpi3RaQ68CUQAqQCjxhjTtq2eRkYCuQBzxljFlZErMqxxcTEEBMTA8DBgwdZsmQJycnJBAQEADB48GD27NlDQkIC999/P9HR0VqrV0o5nYqqwuQCLxhjGgGtgGdEpDEwClhqjAkDlto+Y1v2KPB/gK7AFBFxraBYlZOoV68eQ4cOZcqUKfYEHh0dTW5uLmPGjCEmJobg4GBeffVViyNVSqmyVSHJ3Rhz1Biz0fb+DLADCAa6A9Nsq00DetjedwdmGmMuGGMOAHuBlhURq3Juzz//POvXr+fYsWNMnTqVtm3bUvBoKjc3l+7du/P2229z8OBBiyNVSqmbV+HP3EUkBFgBRAKHjDF+hZadNMZUE5H3gLXGmOm28iTgB2PMN1fbrz5zV6V16NAh7rvvPrZv3w7k3+Lv0aMHAwcOLDKlrlJKVQZXe+ZeoS2LRKQqMAv4H2PM79datZiyIlchIvKEiGwQkQ0ZGRllFaa6Rd1xxx38+uuv7N69mwkTJuDl5cWrr77K3r17Adi3bx+rVq2yT4GrlFKVVYUldxFxJz+xzzDGfGsrPi4iQbblQUC6rTwNqFto8zrAkSv3aYz5yBgTa4yJDQwMLL/g1S0lLCyMESNG8PPPP3PkyBHatWsHwIcffkjbtm2pXbs2jz/+OIsWLSI3N9fiaJVSqqgKSe6S35opCdhhjPlnoUVzgYG29wOBOYXKHxURTxEJBcKA9RURq1KF1apVCze3/E4lo0ePZubMmdxzzz3MnDmTLl26EB4ebp/XXimlKouKGsQmHugPbBORzbayV4DxwFciMhQ4BDwMYIz5VUS+AraT39L+GWOM3gtVlvLx8aF379707t2b8+fPs2DBAo4ePWrvN9+hQwciIiLo3bs38fHxuLpqBw+llDV0EBulysC5c+cYOHAg33//PefPnycoKIhevXrxxBNPEBkZaXV4SiknVSka1CnlrLy9vfn6669JT0/niy++oFWrVnz00Uds3boVgMzMTHbv3m1xlEqpW4Umd6XKUNWqVXn00Uf59ttvSU9PJzExEYBPP/2U8PBw4uLieO+999DeHUqp8qTJXaly4uPjg5eXFwB9+/Zl4sSJXLhwgWHDhlG7dm169Oih3eqUUuVCk7tSFaB27dq88MILbN68ma1btzJ8+HCqVatmb3T3z3/+k9WrV+MsbWCUUtbSBnVKWezs2bPUqVOH06dPExERwZAhQxgwYAA1a9a0OjSlVCWnDeqUqqSqVq1KWloaSUlJ+Pv789JLLxEcHMysWbOsDk0p5aA0uStVCVStWpUhQ4awatUqduzYwfDhw2ndujUAc+fOZdSoUezatcviKJVSjkKTu1KVTEREBBMmTCAoKAiA5ORkJk6cSEREBO3ateM///kPFy5csDhKpVRlpsldqUpu3LhxHD58mPHjx3PkyBH69u1L586drQ5LKVWJaXJXygEEBQUxcuRIdu/ezcKFCxk5ciQA2dnZ9OzZk7lz52q3OqWUnSZ3pRyIi4sLnTt3JiEhAYBdu3axZs0aunfvTmhoKG+88QZHjx61OEqllNU0uSvlwGJiYjh48CCzZs0iIiKCMWPGcMcdd3DgwAGrQ1NKWajEs8KJyG3AHzpLm1KVg7u7O4mJiSQmJrJnzx7mzZtHaGgoAH//+9+pXbs2jzzyCJ6enhZHqpSqKNcdxEZEXIBHgb5AC+AC4AlkAPOBj4wxe8o5zuvSQWyUulxeXh4tW7Zk48aN1KxZkz//+c88+eST1KpVy+rQlFJlpDSD2CwD7gReBmoZY+oaY2oAbYG1wHgR6Vem0SqlSs3V1ZUNGzawcOFCmjdvztixY6lXrx6ff/651aEppcrZjdyW72iMybmy0BhzApgFzBIR9zKPTClVaiJC586d6dy5M7t37+bdd98lLi4OyO8/f+DAAR566CH7GPdKKedw3Zp7cYn9ZtZRSlmrYcOGvPvuuzRs2BCAjz/+mIcffpjw8HCmTJnCuXPnLI5QKVVWbri1vIjEish/RWSjiGwVkW0isrU8g1NKlZ/33nuPb775hoCAAJ555hnq1avH5MmTrQ5LKVUGStIVbgbwKdATeBB4wPZXKeWAXF1d6dmzJ2vWrGHFihXcdddd/P7770B+Y7x9+/ZZHKFS6maVpCtchjFmbrlFopSyhIjQtm1b2rZta59Pfs6cOfTq1YvExERGjBhhf06vlHIMJam5vyoiH4vIn0QkseBVbpEppSqciADQunVrRo0axdKlS2nVqhV33303ixYt4npdZ5VSlcN1+7nbVxSZDkQAvwKXbMXGGDOknGIrEe3nrlTZO3PmDElJSfzjH//Aw8ODXbt24eZW4rGvlFLl5Gr93EuS3LcZY6LKPLIyosldqfJz4cIFUlNTCQ8P5/z58zz44IMMHTqURx55RLvRKWWh0gxiU2CtiDQuw5iUUg7C09OT8PBwAA4dOsTRo0fp06cPERERJCUlcfHiRYsjVEoVVpLk3gbYLCK7tCucUreu8PBwtm3bxqxZs/Dx8eGxxx6jQYMGHDt2zOrQlFI2JXl41rXcolBKORQXFxcSExN56KGHWLhwId999x01a9YEYOXKlcTGxlKlShWLo1Tq1nUjE8eIuc5KN7JOedNn7kpZLysri7p16+Ln58fLL7/M448/jpeXl9VhKeW0SjVxjIgME5E7rtihh4jcKyLTgIFlFahSynH5+/vzww8/EBYWxnPPPUeDBg2YMmUKFy5csDo0pW4pN5LcuwJ5wBcickREtovIfmAP8CdgkjFmajnGqJRyIHfffTfLly9n6dKlhIaG8uyzz+pod0pVsBvuCgdgm/0tADhvjDlVXkHdDL0tr1TlY4whJSWFqKj8XrTDhw+ncePGDBw4EHd3nUxSqdIqi65wGGNyjDFHS5rYReQTEUkXkZRCZdVFZLGI7LH9rVZo2csistfWMr9LSY6llKo8RMSe2C9cuMCaNWt4/PHHCQ8PZ/r06eTl5VkcoVLOqUTJvRSmUrS1/ShgqTEmDFhq+4ytL/2jwP+xbTNFRHSUDKUcnKenJ6tXr+a7777Dz8+P/v37Ex0dzZYtW6wOTSmnUyHJ3RizAjhxRXF3YJrt/TSgR6HymcaYC8aYA8BeoGVFxKmUKl8iwv3338+GDRuYOXMm7u7uBAUFAXD27FmLo1PKeZRkPvd7RSRJRP4hIoNFpLmIeJbi2DWNMUcBbH9r2MqDgcOF1kuzlSmlnISLiwu9e/cmOTmZGjVqYIyhU6dOJCQksHnzZqvDU8rhlaTmPh34DlgL1Af+Sv4kMmVNiikrttWfiDwhIhtEZENGRkY5hKKUKk8Fs9Dl5eWRmJjI2rVradasGX369GHv3r0WR6eU4ypJct9rjPmvMeZrY8wYY0x3Y0yDUhz7uIgEAdj+ptvK04C6hdarAxwpbgfGmI+MMbHGmNjAwMBShKKUspKbmxsjRoxg//79vPLKK8yZM4dGjRqxZMkSq0NTyiFdN7mLyGci8j/AGhF5oQyPPZf/P/jNQGBOofJHRcRTREKBMGB9GR5XKVVJ+fn58eabb7Jv3z5efPFF2rRpA8DWrVvJzs62ODqlHMeN1Nyn2darBfQXkYMiMldEXheRh2/kICLyBbAGCBeRNBEZCowHOonIHqCT7TPGmF+Br4DtwALgGWOM9pdR6hZSq1Yt3nrrLby8vMjJyaFbt26EhYXxySefaPc5pW5AiQaxARARN6Ax0BRoYowZUR6BlZQOYqOU81q9ejUvvPACa9euJSoqiokTJ9K5c2erw1LKcmUyiA2AMSbXGLPVGPN5ZUnsSinn1rp1a1avXs1XX33F2bNn6dKlC8uWLbM6LKUqrYoaxEYppUpFRHj44YfZsWMHn3/+Oe3btwfg+++/58iRYtvcKnXL0uSulHIonp6e9OvXDxHh/PnzDBw4kLCwMMaNG8e5c+esDk+pSkGTu1LKYVWpUoV169aRkJDA2LFjadSoEV999RUlbUuklLPR5K6Ucmh33nknX3/9NcuWLcPPz4/evXuzbds2q8NSylKa3JVSTqF9+/Zs3LiRJUuW0KRJEwD+85//oKNXqluRJnellNNwdXWlQ4cOABw7dowhQ4YQFhbG5MmTycnJsTg6pSqOJnellFOqVasWmzZtIi4ujueff54mTZqwcOFCq8NSqkJocldKOa1GjRqxYMEC5s2bR25uLj169OD48eNWh6VUudPkrpRyaiLCAw88QEpKCkuWLKFmzZoYY/j888/5448/rA5PqXKhyV0pdUvw9PQkPj4egPXr1zNgwAAiIyP54YcfLI5MqbKnyV0pdcuJi4tj8eLFuLq6kpCQQGJiIocOHbI6LKXKjCZ3pdQtqWPHjmzdupW//e1vLFiwgHvuuUdnnFNOQ5O7UuqW5enpycsvv8yOHTtISkrC1dWVnJwcVq5caXVoSpWKJnel1C2vXr169oloPvnkE9q1a0fv3r05evSotYEpdZM0uSulVCEDBw5k3LhxzJkzh4iICN5//30uXbpkdVhKlYgmd6WUKsTLy4u//vWvpKSk0KJFC55++mkee+wxq8NSqkTcrA6gPOXk5JCWlqZ9WZXD8PLyok6dOri7u1sdyi2vQYMGLF68mOnTp9OgQQMAzp49i6urK1WqVLE4OqWuzamTe1paGrfffjshISGIiNXhKHVNxhiysrJIS0sjNDTU6nAU+QPg9O/f3/551KhRLFiwgA8++ICOHTtaGJlS1+bUt+X/+OMP/P39NbErhyAi+Pv7652mSqxnz564uLjQqVMn+vfvrzPOqUrLqZM7oIldORT9vVZu99xzD1u3bmX06NF8+eWXREREsGjRIqvDUqoIp0/uVnvnnXdo1KgRffv2tToUxo4dy8SJE0u1j+XLl/PAAw9cc51Tp04xZcoU++cjR47Qq1evUh23smrduvVNbVcW/y2UNby8vHj99dfZvHkzzZs358477wTyH6soVVloci9nU6ZMYf78+cyYMeOG1s/NzS2T4xpjLOu+c2Vyr127Nt98840lsZSXgpHMVq9ebXEkyiqNGzdm0aJF3HnnnRhjSExMZOTIkWRnZ1sdmlKa3MvTU089xf79++nWrRuTJk3ixIkT9OjRgyZNmtCqVSu2bt0K5NfinnjiCTp37syAAQNISEiwL2vWrBmvvfYaAGPGjOHjjz/m7NmzdOjQgZiYGKKiopgzZw4AqampNGrUiKeffpqYmBgOHz7Mm2++SXh4OB07dmTXrl3Fxvn1118TGRlJ06ZNadeuHZDfXmHw4MFERUXRrFkzli1bVmS7K2ufkZGRpKamMmrUKPbt20d0dDQjRowgNTWVyMjIa+536tSpJCYm0rVrV8LCwnjppZeue37bt2/PyJEjadmyJQ0bNrSPKjZ16lSeffZZ+3oPPPAAy5cvB6Bq1aqMHDmS5s2b07FjR9avX0/79u2pX78+c+fOBfIT94gRI2jRogVNmjThww8/BPLvWtxzzz306dOHqKgo+/4KTJgwgaioKJo2bcqoUaMA+Pe//02LFi1o2rQpPXv25Ny5c9f9Xsrx5OTkUL16dSZMmEBISAhxcXF07drV/ps6ceIEkyZNYtq0acydO5dVq1bx66+/cubMGYsjV87KqVvLX6lgBKrCHnnkEZ5++mnOnTtHQkJCkeWDBg1i0KBBZGZmFrm1XJAwruaDDz5gwYIFLFu2jICAAIYNG0azZs2YPXs2P/74IwMGDGDz5s0AJCcns2rVKqpUqcL48eNZuXIlISEhuLm58fPPPwOwatUq+vXrh5eXF//973/x8fEhMzOTVq1a0a1bNwB27drFp59+ypQpU0hOTmbmzJls2rSJ3NxcYmJiaN68eZE4X3vtNRYuXEhwcDCnTp0C4F//+hcA27ZtY+fOnXTu3Jndu3df8/sWGD9+PCkpKfbvlpqaal92rf1u3ryZTZs24enpSXh4OMOGDaNu3brXPFZubi7r169n/vz5jBs3jiVLllxz/ezsbNq3b8/f//53HnroIUaPHs3ixYvZvn07AwcOpFu3biQlJeHr68svv/zChQsXiI+Pp3PnzkD+bGIpKSlFWrP/8MMPzJ49m3Xr1uHt7c2JEycASExM5PHHHwdg9OjRJCUlMWzYsBs6j8pxeHh4kJSUxJAhQ3j//ffJyMjgxIkT9saRBw4cYPjw4UW2++yzz+jfvz9r166lV69eVKtWjerVq9v//uUvf6Fp06akpaWxatWqy5ZVr14dX19fXFy0jqaKuqWSu9VWrVrFrFmzALj33nvJysri9OnTAHTr1s3ed7Zt27a88847hIaGcv/997N48WLOnTtHamoq4eHh5OTk8Morr7BixQpcXFz47bffOH78OJA/jGarVq0AWLlyJQ899BDe3t72YxQnPj6eQYMG8cgjj5CYmGiPtSAJRUREUK9evRtO7tc7B1fbb4cOHfD19QXyb3kePHjwusm9IN7mzZtfdhFxNR4eHnTt2hWAqKgoPD09cXd3Jyoqyr79okWL2Lp1q/1RwunTp9mzZw8eHh60bNmy2G5qS5YsYfDgwfZzXb16dQBSUlIYPXo0p06d4uzZs3Tp0uW6MSrHFR8fb59WtrBmzZpx4sQJTp48ednfuLg4AHx8fOjSpYu9/MCBAyQnJzNgwAAA1qxZw5/+9Kci+125ciVt2rRhzpw5vPHGG1SvXp2AgAD8/f0JCAjgqaeeokaNGhw5coT09HT8/f3x9/e3/06V87qlkvu1atre3t7XXB4QEHDdmvr1FNfgpqB19G233WYva9GiBRs2bKB+/fp06tSJzMxM/v3vf9tr3TNmzCAjI4Pk5GTc3d0JCQmx1xAK76fw/q/lgw8+YN26dXz//fdER0ezefPmG2oc5Obmdtlz/RvpwnWt/Xp6etrfu7q63lD7g4JtCq9/rbjc3d3t58TFxcW+vYuLi317YwzvvvtukUS8fPnyIue38Pcq7lwPGjSI2bNn07RpU6ZOnVrq35ByTC4uLlSrVo1q1apRv379IssbN25MUlLSVbdPSEhg+/btl10YnDhxgrCwMCC/kV9AQAAnTpxg3759ZGZmcvr0afr160eNGjWYPn06I0eOtO+vSpUq+Pv7k5ycTI0aNfj2229ZunSp/aKg4CKgU6dO9sl03NzctDeHA7mlkrvV2rVrx4wZMxgzZgzLly8nICAAHx+fIut5eHhQt25dvvrqK8aMGUNGRgYvvvgiL774IpBfk6xRowbu7u4sW7aMgwcPXvV4gwYNYtSoUeTm5jJv3jyefPLJIuvt27ePuLg44uLimDdvHocPH7bHeu+997J7924OHTpEeHg4a9assW8XEhLCd999B8DGjRs5cOAAALfffvtVnyVebb8bN2686nkbMGAAzz77LC1btrzqOoWFhIQwZcoULl26xG+//cb69etvaLsCXbp04f333+fee+/F3d2d3bt3ExwcfM1tOnfuzGuvvUafPn3st+WrV6/OmTNnCAoKIicnhxkzZlx3P0oV57bbbqNRo0ZXXd6lS5ciF6M5OTm4uroC0KtXL8LCwsjMzCQrK4usrCwyMzPt///Zvn07M2fO5OTJk/YLcBEhJycHgOeee45PPvnksuRfu3Zte0PhRYsWkZ6eTmBgIDVq1CAwMJDAwMDLLthVxdLkXoHGjh3L4MGDadKkCd7e3kybNu2q67Zt25alS5fi7e1N27ZtSUtLo23btgD07duXBx98kNjYWKKjo4mIiCh2HzExMfTu3Zvo6Gjq1atn3/5KI0aMYM+ePRhj6NChA02bNiUiIoKnnnqKqKgo3NzcmDp1apF/qD179uSzzz4jOjqaFi1a0LBhQwD8/f2Jj48nMjKS++67j2eeeca+zdNPP33d/V5p69atBAUFXXOdwuLj4wkNDSUqKorIyEhiYmJueFuAxx57jNTUVGJiYjDGEBgYyOzZs6+5TdeuXdm8eTOxsbF4eHiQkJDA3/72N15//XXi4uKoV68eUVFR2oBKVZjCQxjXr1+/2DsGBUaPHs3o0aPJy8vj5MmTZGVlcfLkSfvFwX333Yevr+9lFwdpaWn27d977z3mzZt32T5DQkLsF/zDhw9n37599qRfo0YNGjRowIMPPghAZmYmt99+u14MlCFxlr6ZsbGxZsOGDZeV7dix45pXu6ry+/333xk6dChff/211aFUGP3dKkdz+vRpjh8/Tnp6OhkZGWRkZODq6srQoUOB/Iv6n3/+2b4sNzeX1q1b2xsLR0VFkZKSgo+Pj/0CoEOHDrzxxhsAfPrpp7i7u+udgWKISLIxJvbKcq25q0rNx8fnlkrsSjkiX19ffH197XfvrlR43AtjDKdOneL8+fP2spEjR5KammpP/unp6Vy8eNG+fPjw4faePAX69u3L9OnTAejUqRM+Pj7UrFnT/mrRogXNmzfHGMP58+dvuUaElTq5i0hX4G3AFfjYGDPe4pCUUkqVgojYGxcW6Nev3zW32b9/vz3xFyT/gl4rOTk55OXlsXPnTn766SeysrIAeOmll2jevDm///47fn5+VK1a9bLkP2jQILp168a5c+dYuHDhZcsKj1/hqCptchcRV+BfQCcgDfhFROYaY7ZbG5lSSqmKVHAxUNydAXd3d3788Uf754sXL5KRkWFvc+Di4sJbb73F8ePH7a/du3fbJ/05cOCAvUttAW9vbz788EP69evHgQMHmDBhwmXJv2bNmkRGRtq77lZGlTa5Ay2BvcaY/QAiMhPoDpQouV+ti5JSlZGztIFRyioeHh6X9Uq5/fbb7SNGFufOO+8kOTn5suR//PhxwsPDgfy5Mb755huysrIu+/c5d+5cHnzwQZYsWcIzzzxDrVq1qFWrFkFBQdSqVYv+/fsTHBzMmTNn7DOUVuSAQ5U5uQcDhwt9TgPiSrIDLy8vsrKydNpX5RAK5nP38vKyOhSlbhleXl7X7FETHx9vbwSYkZHBsWPHOH78uH3ckapVqxIdHc2xY8fYtGkT8+fPtw9YFRwczBdffMGTTz6Jm5sbkyZNumxo7PJUmZN7cdn4smqNiDwBPAFwxx13FFm5Tp06pKWl6ZzLymF4eXlRp04dq8NQSl3Bzc2NoKCgIt1yW7VqxZdffnlZ2dmzZ+0X6W3atOGdd97h6NGjNGvWrOLirbAjlVwaUHjs0TrAkcIrGGM+Aj6C/K5wV+7A3d292KFClVJKqfJSuEFe48aNady4cYXHUJlnHPgFCBORUBHxAB4F5lock1JKKVXpVdqauzEmV0SeBRaS3xXuE2PMrxaHpZRSSlV6lTa5Axhj5gPzrY5DKaWUciROM/ysiGQAxc+gcvMCgMwy3uetRs9h6ek5LD09h6Wn57D0yuMc1jPGBF5Z6DTJvTyIyIbixuxVN07PYenpOSw9PYelp+ew9CryHFbmBnVKKaWUugma3JVSSikno8n92j6yOgAnoOew9PQclp6ew9LTc1h6FXYO9Zm7Ukop5WS05q6UUko5GU3uxRCRriKyS0T2isjVpxNSxRKRuiKyTER2iMivIvIXq2NyVCLiKiKbROQ7q2NxVCLiJyLfiMhO22/yLqtjcjQi8rzt33KKiHwhIjq70XWIyCciki4iKYXKqovIYhHZY/tb7Vr7KA1N7lcoNI/8fUBj4E8iUvEDAzu2XOAFY0wjoBXwjJ7Dm/YXYIfVQTi4t4EFxpgIoCl6PktERIKB54BYY0wk+SOGPmptVA5hKtD1irJRwFJjTBiw1Pa5XGhyL8o+j7wx5iJQMI+8ukHGmKPGmI2292fI/59p8LW3UlcSkTrA/cDHVsfiqETEB2gHJAEYYy4aY05ZGpRjcgOqiIgb4M0Vk3ipoowxK4ATVxR3B6bZ3k8DepTX8TW5F1XcPPKamG6SiIQAzYB1FofiiCYDLwGXLI7DkdUHMoBPbY83PhaR26wOypEYY34DJgKHgKPAaWPMImujclg1jTFHIb8SBNQorwNpci/quvPIqxsjIlWBWcD/GGN+tzoeRyIiDwDpxphkq2NxcG5ADPC+MaYZkE053gp1Rrbnwt2BUKA2cJuI9LM2KnU9mtyLuu488ur6RMSd/MQ+wxjzrdXxOKB4oJuIpJL/aOheEZlubUgOKQ1IM8YU3Dn6hvxkr25cR+CAMSbDGJMDfAu0tjgmR3VcRIIAbH/Ty+tAmtyL0nnkS0lEhPxnnDuMMf+0Oh5HZIx52RhTxxgTQv5v8EdjjNaWSsgYcww4LCLhtqIOwHYLQ3JEh4BWIuJt+7fdAW2UeLPmAgNt7wcCc8rrQJV6ylcr6DzyZSIe6A9sE5HNtrJXbFP4KlXRhgEzbBfr+4HBFsfjUIwx60TkG2Aj+T1hNqGj1V2XiHwBtAcCRCQNeBUYD3wlIkPJv2h6uNyOryPUKaWUUs5Fb8srpZRSTkaTu1JKKeVkNLkrpZRSTkaTu1JKKeVkNLkrpZRSTkaTu1JKKeVkNLkrpZRSTkaTu1LqholIFRH5yTbPfEjhuapLuB8PEVlhm2VMKVXGNLkrpUpiCPCtMSavNDuxTae8FOhdJlEppS6jyV0pBYCILBORTrb3b4jIO8Ws1pdixsMWkfq2KVVb2Gr0O23Tq6aIyAwR6SgiP4vIHhFpadtstm1/SqkypsldKVXgVeD/ikhfoBnwfOGFtrHZ6xtjUq8oDyd/BsDBxphfbMUNgLeBJkAE0AdoA7wIvGJbJwVoUS7fRKlbnD7vUkoBYIxZYZv1azjQvphb7wHAqSvKAsmvyfe8YoKlA8aYbQAi8iuw1BhjRGQbEGI7Xp6IXBSR240xZ8r+Gyl169Kau1IKABGJAoKAC1dJtucBryvKTgOHyZ8JsLALhd5fKvT5EpdXKjyBP242ZqVU8TS5K6UQkSBgBtAdyBaRLleuY4w5CbiKSOEEfxHoAQwQkT4lPKY/kGGMybnpwJVSxdLkrtQtTkS8gW+BF4wxO4DXgbFXWX0R+c/O7Ywx2cADwPMi0r0Eh74HmF/igJVS16XzuSulbpiINAOGG2P6l8G+vgVeNsbsKn1kSqnCtOaulLphxphNwDIRcS3Nfmwt72drYleqfGjNXSmllHIyWnNXSimlnIwmd6WUUsrJaHJXSimlnIwmd6WUUsrJaHJXSimlnIwmd6WUUsrJaHJXSimlnMz/A3a9l6b0iSefAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import flopy\n",
    "\n",
    "nrow, ncol = 1, int(L)\n",
    "M0 = flopy.modflow.Modflow(modelname='model', exe_name='../mf2005.exe')\n",
    "flopy.modflow.ModflowDis(M0, nlay=1, nrow=nrow, ncol=ncol, nper=1, delr=1., delc=1., top=0., botm=-b, steady=True, \n",
    "                         perlen=1., nstp=1)\n",
    "flopy.modflow.ModflowBas(M0, ibound=np.hstack([np.ones([nrow, ncol-1], dtype=int), -1*np.ones([1,1])]), \n",
    "                         strt=BC1h*np.ones([nrow, ncol], dtype=float))\n",
    "flopy.modflow.ModflowLpf(M0, hk=K, vka=-999., ss=-999., sy=-999., ipakcb=53)\n",
    "flopy.modflow.ModflowRch(M0, nrchop=1, rech=R, ipakcb=53)\n",
    "flopy.modflow.ModflowPcg(M0, hclose=1e-6, rclose=1e-6)\n",
    "flopy.modflow.ModflowOc(M0, stress_period_data={(0,0): ['save head', 'save budget']})\n",
    "M0.write_input()\n",
    "success, buff = M0.run_model(silent=True)\n",
    "H0 = flopy.utils.binaryfile.HeadFile('model.hds').get_data()[0,0,:]\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "plt.subplots(figsize=[8,2])\n",
    "plt.plot(X/1000., H0, 'k--', mfc='none', label='forward solution, numerical')\n",
    "plt.xlabel('$x$ (km)')\n",
    "plt.ylabel('$h$ (m)')\n",
    "plt.legend(loc=3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgAAAACdCAYAAADCH8BOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoNElEQVR4nO3de5yOdf748dd7Do4ROYXBGBFjzAwmp3EsyqkcliIVWllKqTbnX6ms3fqykpJDJ5EIq2FRbTFtabcyNnIWzTjErjHEMDSn9++P+3Y3zAwzmLnmnvv9fDw83Pfncx3e9zX3zPW+Ptfn+nxEVTHGGGOMb/FzOgBjjDHGFD5LAIwxxhgfZAmAMcYY44MsATDGGGN8kCUAxhhjjA+yBMAYY4zxQQFOB1CYKleurMHBwU6HYYwxxhSKzZs3H1fVKjnVeXUCICJdgVcBf+AtVX3pcssHBwcTFxdXKLEZY4wxThORA7nVee0tABHxB2YD3YBQYKCIhDoblTHGGOMdvDYBAFoA+1T1J1VNBZYCvRyOyRhjjPEK3nwLoCZwKMv7w0BLh2Ixxpg8SUpKYuDAgdnKhw8fTr9+/Th8+DAPP/xwtvrRo0fTo0cP9u7dy6hRo7LVjx8/nttvv52tW7cyZsyYbPUvvvgirVq14ptvvuG5557LVj9t2jQiIiLYsGEDL72U/W7q66+/ToMGDVi7di2vvvpqtvp33nmHoKAgVqxYwfz587PVL1myhEqVKrFo0SIWLVqUrT4mJoYyZcowb948/va3v2Wr/8c//gHAzJkzWbdu3UV1pUuXZtWqVQD85S9/ITY29qL6m266iaVLlwLw7LPP8u23315UX7NmTd59910AnnnmGX744YeL6uvXr8/s2bMBeOyxx/jxxx8vqg8PD2f69OkADB06lJ9//vmi+pYtWzJlyhQABgwYwIkTJy6q79SpExMmTMj2mQuaNycAkkNZtokNRGQ4MBygdu3a2VZIS0vj8OHDnD9//roHaExxV6pUKYKCgggMDHQ6lCJNVdm3bx/169cnMzOTM2fOZFsmLS0NINf69PT0PNVnZGRctj49PT3H+oyMjMvWZ2Zm5qk+LS3tsvWpqak51l+Yl+bXX3/Nsf6C8+fPZ6u/sO3c6kuVKnXZ+pSUFM/rc+fOXbY+JSUlW/25c+cuW5/1/HKl+sIk3joZkIi0Bp5X1bvc7ycAqOpfclsnKipKL+0EGB8fT7ly5ahUqRIiOeUUxpicqCpJSUkkJydTt25dp8MpspKSkhg2bBjr169n586dBAUFOR2S8SEisllVo3Kq8+Y+AJuA+iJSV0RKAAOA1fndyPnz5+3kb8xVEBEqVapkrWeXERsbS3h4OGvXruWFF16gRo0aTodkjIfXJgCqmg6MAj4FdgHLVHXH1WzLTv7GXB373cmZqjJu3DjuuOMOypUrx7fffstTTz2Fn5/X/sk1xZBXfxtVdZ2qNlDVeqo61el4nBITE8POnTvzvd4XX3xBz54987XOzJkzL7of1r17d3755Zd877swzZ07l4ULFwKwYMECjhw54qkbNmzYFY9dmzZtAEhISOCDDz4ouEALQdZjkV833HDDdY6m+BIRkpKSeOSRR9i8eTNNmzZ1OiRjsvHmToAGV4ecmJgYevbsSWho3odBuNAhKL9mzpzJAw88QJkyZQCy9cYtikaMGOF5vWDBAsLCwjxNsW+99dYV1//Xv/4F/JYA3H///QUTaAFLT0+/6FiY60tVeeutt4iKiqJp06bMmzcPf39/p8MyJlde3QJQHCQkJNCwYUMGDx5MeHg4/fr181xhb968mQ4dOtC8eXPuuusujh49CkDHjh2ZOHEiHTp04OWXX2b16tWMGTOGyMhI9u/fT8eOHT0jHh4/fpwLwx8vWLCA/v37c/fdd3PnnXcCcPr0afr06UNoaCgjRozw9KYdOXIkUVFRNG7cmMmTJwMwa9Ysjhw5QqdOnejUqRPgGl3x+PHjAMyYMYOwsDDCwsKYOXOm5/M1atSIRx55hMaNG3PnnXde1GP2guXLlxMWFkZERATt27cHXL2Sx4wZw2233UZ4eDjz5s0DXC0XHTt2pF+/fjRs2JBBgwZ5ehCPHz+e0NBQwsPDeeaZZwB4/vnnmT59OitWrCAuLo5BgwYRGRnJuXPnPMdqzpw5jB071hPPggULePzxx4HfrnzHjx/PV199RWRkJK+88grt2rVjy5YtnnWio6OzPT6U1eXiznoc4+Li6Nixoyf2wYMHc+eddxIcHMzKlSsZO3YsTZo0oWvXrp6e43n5rrz66queYwGwb98+OnfuTEREBM2aNWP//v2cOXOGO+64g2bNmtGkSRPPo1Xm8o4dO0bfvn0ZPny45xE4O/mbIk9VfeZf8+bN9VI7d+7MVlaY4uPjFdCNGzeqqurQoUN12rRpmpqaqq1bt9Zjx46pqurSpUt16NChqqraoUMHHTlypGcbgwcP1uXLl3ved+jQQTdt2qSqqomJiVqnTh1VVX333Xe1Zs2ampSUpKqqsbGxWrJkSd2/f7+mp6dr586dPdu5sEx6erp26NBBt27dqqqqderU0cTERM++LryPi4vTsLAwPXPmjCYnJ2toaKj+5z//0fj4ePX399fvv/9eVVX79++vixYtynYcwsLC9PDhw6qqevLkSVVVnTdvnk6ZMkVVVc+fP6/NmzfXn376SWNjY7V8+fJ66NAhzcjI0FatWulXX32lSUlJ2qBBA83MzLxoO5MnT9Zp06ZlOzZZ3x87dkzr1avnKe/atat+9dVXqqpatmxZz/Hq0aOHZ5kFCxbo6NGjVVV1z549mtP3K6vc4r70uG7atEk7dOjgiT06OlpTU1N1y5YtWrp0aV23bp2qqvbu3Vs/+uijfH1Xsh6LFi1a6MqVK1VV9dy5c3r27FlNS0vTU6dOqarru1OvXj3P8bxwHC7l9O+Q09asWaNVq1bVkiVL6vTp0zUjI8PpkIzxAOI0l3Oi3QK4xIUrr6zuvfdeHn30UVJSUujevXu2+iFDhjBkyBCOHz9Ov379Lqr74osvrrjPWrVqER0dDcADDzzArFmz6Nq1K9u3b6dLly6A62q4evXqnnXuu+++fHyq33Tp0oWbbrrJ875FixaEhIQAMHDgQDZu3Ei/fv1YtmwZ8+fPJz09naNHj7Jz507Cw8Nz3e7GjRvp06cPZcuWBaBv37589dVX3HPPPdStW5fIyEgAmjdvTkJCQrb1o6OjGTJkCPfeey99+/YFXAN//PDDD6xYsQKAU6dO8eOPP1KiRAlatGjheZwqMjKShIQEWrVqRalSpRg2bBg9evTIV/+GKlWqEBISwjfffEP9+vXZs2eP52eSm/79+zNlyhSmTZvGO++8w5AhQ664n5zibtu27WXX6datG4GBgTRp0oSMjAy6du0KQJMmTUhISGDPnj35/q4kJyfz888/06dPH+C356TT0tKYOHEiX375JX5+fvz888/873//4+abb77iZ/NFq1evplevXjRp0oT169cTFhbmdEjG5JklAEXApT2pRQRVpXHjxvz73//OcZ0LJ9qcBAQEeJryL31E69L1ctp3fHw806dPZ9OmTVSsWJEhQ4Zc8VEvvcx4EiVLlvS89vf3z/EWwNy5c/n2229Zu3YtkZGRbNmyBVXltdde46677rpo2S+++CLbNtPT0wkICOC7775j/fr1LF26lNdff50NGzZcNu6s7rvvPpYtW0bDhg3p06fPFXu4lylThi5durBq1SqWLVuWp4mmcoobLv8zu7COn58fgYGBnrj8/PxIT0+/qu9Kbj+vxYsXk5iYyObNmwkMDCQ4ONge88tBSkoKZcqUoVu3bsyYMYNHH330op+tMd7A+gBc4osvvsj279FHHwVcf/Bzqr9w5Ve5cuVsdXlx8OBBzx/vJUuW0LZtW2699VYSExM95WlpaezYkfNTjuXKlSM5OdnzPjg4mM2bNwN4rp5z89133xEfH09mZiYffvghbdu25fTp05QtW5Ybb7yR//3vf3z88ce57uuC9u3bExMTQ0pKCmfPnuWjjz6iXbt2efr8APv376dly5a8+OKLVK5cmUOHDnHXXXcxZ84cz33uvXv3cvbs2Vy3cebMGU6dOkX37t2ZOXPmRffnrxQ/uFotYmJiWLJkSY5XzTmtO2zYMJ544gluu+02T8vKd999x0MPPZTXjw5c/DPLaRjUy8nPd+WC8uXLExQURExMDOAafS0lJYVTp05RtWpVAgMDiY2N5cCBXCcS80mpqak899xzNGrUiBMnThAYGMhTTz1lJ3/jlSwBKAIaNWrEe++9R3h4OCdOnGDkyJGUKFGCFStWMG7cOCIiIoiMjPT0Rr/UgAEDmDZtGk2bNmX//v0888wzzJkzhzZt2ng6luWmdevWjB8/nrCwMOrWrUufPn2IiIigadOmNG7cmIcffviipvDhw4fTrVs3TyfAC5o1a8aQIUNo0aIFLVu2ZNiwYfl69GnMmDE0adKEsLAw2rdvT0REBMOGDSM0NJRmzZoRFhbGH/7wh8s+vZCcnEzPnj0JDw+nQ4cOvPLKK9mWGTJkCCNGjPB0AsyqYsWKhIaGcuDAAVq0aJFt3fDwcAICAoiIiPBsu3nz5pQvX56hQ4d6ljt48CClS5fO82cHmDx5MqNHj6Zdu3b57jyWn+9KVosWLWLWrFmEh4fTpk0b/vvf/zJo0CDi4uKIiopi8eLFNGzYMF+xFGc7duygVatWTJkyhY4dOxIQYA2oxrt57VDAVyOnoYB37dpFo0aNHIrI1Uu+Z8+ebN++3bEYzNU7cuQIHTt2ZPfu3Z5BXsaMGcODDz542T4TxYnTv0MFLTMzk7/+9a88++yzlC9fnvnz59O7d2+nwzImT4rrUMDGOGrhwoW0bNmSqVOnXjTC27Rp03zm5O8LRIQNGzbQvXt3tm3bZid/U2xYC0Axv3oxpqAVx9+hzMxM3njjDXr27ElwcDApKSmULl3ahj42XsdaAIwxJo9++uknOnXqxOOPP87bb78NuDoA28nfFDfWiwXXI1H2y21M/hWnFsTMzEzmzp3L2LFj8ff3591332Xw4MFOh2VMgfH5FoBSpUqRlJRUrP6QGVMYVJWkpCTPIELebvr06Tz22GO0bt2abdu2MWTIELswMMVakW4BEJHngUeARHfRRFVd566bAPweyACeUNVPr2YfQUFBHD58mMTExCsvbIy5SKlSpTwjG3qjjIwMjh49SlBQECNGjKBWrVoMGDDATvzGJxTpBMDtFVWdnrVAREKBAUBjoAbwuYg0UNWM/G48MDCQunXrXp9IjTFeY9euXfz+97/n9OnTfP/995QvX56BAwc6HZYxhcZbbwH0Apaq6q+qGg/sA7KP3GKMMZdITU1lypQpREREsHfvXsaPH2+D+hif5A3f+lEi8hAQB/xRVU8CNYFvsixz2F1mjDG5Onz4MN26dWP79u3ce++9zJo1i2rVqjkdljGOcLwFQEQ+F5HtOfzrBcwB6gGRwFHgrxdWy2FTOfbiE5HhIhInInF2n98Y33Shk+/NN99MSEgIq1ev5sMPP7STv/FpjrcAqGrnvCwnIm8Ca9xvDwO1slQHAUdy2f58YD64BgK6+kiNMd7ok08+4YUXXmDNmjVUqlSJVatWOR2SMUWC4y0AlyMi1bO87QNcGDB/NTBAREqKSF2gPvBdYcdnjCm6jh07xv3330+3bt345ZdfOHr0qNMhGVOkON4CcAX/JyKRuJr3E4A/AKjqDhFZBuwE0oHHruYJAGNM8bRgwQL++Mc/kpyczOTJk5kwYYJN2WvMJYp0AqCqD16mbiowtRDDMcZ4iXXr1hEWFsacOXMIDQ11OhxjiqQifQvAGGPyIiUlhQkTJrBnzx4AZs+eTWxsrJ38jbmMIt0CYIwxV7Jq1Sqeeuop4uPjqVatGrfeeitVqlRxOixjijxrATDGeKWEhAR69+5N7969KV26NOvXr+fJJ590OixjvIYlAMYYrzRnzhw+++wzXnrpJbZs2cLtt9/udEjGeBXxpVnwoqKiNC4uzukwjDFX6fPPP6dUqVK0bduWs2fPcvz4cerUqeN0WMYUWSKyWVWjcqqzFgBjTJGXkJDAvffeS5cuXfjTn/4EQNmyZe3kb8w1sATAGFNknTt3jhdeeIFGjRqxZs0aXnjhBWJiYpwOy5hiwZ4CMMYUWcuXL+f555+nX79+zJgxg1q1al15JWNMnlgCYIwpUnbv3s1PP/1E9+7deeCBB2jQoAGtWrVyOixjip0r3gIQEf/CCMQY49t++eUXnn76aZo0acLIkSNJT0/Hz8/PTv7GFJC89AF4q8CjMMb4rPT0dObOnUv9+vWZOXMmQ4cOZdOmTQQEWAOlMQUpLwmAv4i8eGmhiJQTkeUFEJMxxods2rSJkSNHEhoayubNm5k/fz5Vq1Z1Oixjir28JAAPA61FZNiFAhGJAOKAXQUVmDGm+NqxYwdvvvkmAK1bt2bjxo188cUXNG3a1OHIjPEdV0wAVDUd6AuMFJHuIvIIsArXFLzPFXSAxhgvtngxBAeDnx8EB3NqzhxGjRpFREQEEydOJDk5GYDo6GhExNlYjfExeekE+AYwCPg/4F3gPqClqn5+PQIQkf4iskNEMkUk6pK6CSKyT0T2iMhdWcqbi8g2d90ssb8cxhQ9ixfD8OFw4ACowoEDBD76KL+88QbDhg1j9+7dlCtXzukojfFZeellswUIB5oAgUAD4E0R2QZsU9Wl1xjDdlwtDPOyFopIKDAAaAzUAD4XkQaqmgHMAYYD3wDrgK7Ax9cYhzHmepo0CVJSLioqA7xbvTqBc+c6E5MxxuOKCYCqzs/6XkSC+C0h6AFcUwKgqrvc2720qhewVFV/BeJFZB/QQkQSgPKq+m/3eguB3lgCYEyRoapw8CA5Nc0FHj1a6PEYY7LL91DAqnpYVdep6suq+mBBBOVWEziU5f1hd1lN9+tLy3MkIsNFJE5E4hITEwskUGPMb9577z3CwsI4kNtEY7VrF25Axpgc5elBWxEpCYwFWgNngf8Ay1R1fx7X/xy4OYeqSaq6KrfVcijTy5TnyN2CMR9cswFeIVRjzFVKTU0lOjqauLg4/Pz82DVoEHU++gjJehugTBmYOtW5II0xHnltAZgG1AEm4LrirgS8KyIP5GVlVe2sqmE5/Mvt5I97P1kH/g4CjrjLg3IoN8Y44Pjx4wCUKFGCMmXKEB0dzbFjx+j2/vvI/PlQpw6IuP6fPx8GDXI4YmMM5D0BiACeUdWtQAtVfQboAjxWYJHBamCAiJQUkbpAfeA7VT0KJItIK3fv/4dwPZZojClEBw8e5PHHH6dWrVps3boVgH/+859s3LiRSpUquRYaNAgSEiAz0/W/nfyNKTLyOtbmNH5LFo6LyAxgN5B5rQGISB/gNaAKsFZEtqjqXaq6Q0SWATuBdFzjDmS4VxsJLABK4+r8Zx0AjSkkJ0+eZMaMGfzpT38CYODAgdx0000OR2WMyS/R3DrqXFhApDXwjboXFJESQD9c9/SXqqrXNL9HRUVpXFyc02EY47VUlaZNm7J161YaNmzIe++9R4sWLZwOyxiTCxHZrKpROdXlpQVgMDBbRPYCnwCfqOoH1zNAY0zRdfbsWfr27cvcuXOpW7cub7zxBmXKlCEyMtLp0Iwx1yAv4wCMABCRhkA3YIGI3AjE4koIvs7SNG+MKSbOnj3L7NmzGTduHADffvstdevWpU2bNg5HZoy5HvI836aq7sZ13/8VESkNdAL6AzOAHJsXjDHeR1V55513GD9+PMePH6dixYrExMTQvn17p0MzxlxHeU4AROR2XHMC/IJr+N4fcD0Z8GvBhGaMKUzJycmUKlWKgIAAvv/+eyIjI5k0aRIdO3Z0OjRjTAHIz0iA7wNrcI2/HwI8B+woiKCMMYXn5MmTvPTSS9SqVYuXX34ZEWHq1Kl89tlndvI3phjLcwsAsE9VP3K/Xl4QwRhjCs+RI0eYOXMm8+bN4/Tp03Tr1o3OnTsDcOONNzocnTGmoF0xAXBPtvMf4N8i8kdV/WvBh2WMKWhPPvkkf/vb37j77ruZMGECLVu2dDokY0whysstgPfcy90MPCgiB0RktYhMEZH+BRueMeZ6iYuLQ0R44YUXAJg1axbbt28nJibGTv7G+KC8PAa4Hlh/4b2IBAChuIYHbondDjCmyFJVNmzYwMCBA7kwG+bBgwcBuPnmm7n55pzm6DLG+IJrnQ1wUQHGZoy5Ru+//z4PPfQQAL179+btt9+2YXuNMcC1zwZoM3sYU4QkJyczfvx4PvzwQwCio6OZOXMmZ86c4aOPPrKTvzHG41pnAxxVYJEZY/LsyJEj/O53v+PGG2/k5ZdfZsOGDQCEhIQwevRoypYt63CExpiixvHZAI0x12bSpEn8+c9/BuCOO+6gb9++jBw50uGojDFFXZ4SAFVdk+Vtf36bDfCanwJwP0nwPNAIV+tCnLs8GNgF7HEv+k2WeQma89t0wOuA0XqlaQ2NKSbS09OZMWMGJUuWZPTo0XTt2pUzZ84wYsQIGjVq5HR4xhgvkZ+BgABQ1VTges4GuB3oC8zLoW6/qkbmUD4HGI5rVMJ1QFfg4+sYkzFFzsmTJ3nzzTeZPXs2Bw8epH///jzxxBO0a9eOdu3aOR2eMcbL5Gco4AKhqrtUdc+Vl3QRkepAeVX9t/uqfyHQu6DiM6YoWLBgATVq1GDcuHHUrFmThQsXsnjxYkTE6dCMMV4q3y0AhayuiHwPnAb+n6p+BdTE9STCBYfdZcYUG6dPn2bJkiUEBQXRo0cP6tWrx6BBgxg1ahSRkZFOh2eMKQYKJQEQkc9x9Rm41CRVXZXLakeB2qqa5L7nHyMijYGcLnlyvf8vIsNx3S6gdu3a+QvcmEL29ddfs3z5cl599VXANVxvjx49rJnfGHPdFUoCoKqdr2KdX4Ff3a83i8h+oAGuK/6gLIsGAUcus535wHyAqKgo6yhoiqwPPviAQYN+G1rjxRdfZNKkSQ5GZIwpzhzvA5AbEakiIv7u1yFAfeAnVT0KJItIK3HdAH0IyK0VwZgia8+ePdxzzz1MnDgRgDZt2tCvXz8SEhJQVZ599ln8/Irsr6gxxss53gdARPoArwFVgLUiskVV7wLaAy+KSDqQAYxQ1RPu1Uby22OAH2NPABgvkZycTExMDK+99hqbNm0CoGTJkgAEBwezfLlNrWGMKRziS4/PR0VFaVxcnNNhGB+jqqgqfn5+/PnPf2bSpEnUq1eP4cOHM3DgQGrVquV0iMaYYkpENqtqVE511r5oTAGJj48nJCSEkJAQ5s1zDXNx++23s2HDBvbu3cvYsWPt5G+McYwlAMZcZx9++CFVq1YlJCSE+Ph4atSoQalSpQBo1aoVnTp1snv7xhjH2V8hY65RWloaCxYs4Pz58wDs3r2blJQURo8eTWxsLF9//TVDhw51OEpjjLmY9QEw5ipkZmayfv16BgwYwIkTrr6pU6dOZeLEiaSnp+Pn52dX+cYYx12uD4DjTwEY440GDRrE0qVLAWjevDmjRo1iwIABAAQE2K+VMabos0sUYy4jMzOThQsX0rBhQ8qWLUtSUhIA9913HwsWLODEiRPExcUxZMgQz31+Y4zxBnapYkwOMjMzGTp0KAsXLvSUderUiZMnT1KpUiV69+7tXHDGGHMdWAJgDJCamsrGjRuJjY1l6NCh1K5dm4ULFxIZGUmXLl148sknqVGjhtNhGmPMdWMJgPFZp06d4vPPP2fNmjWsXr3a05mvefPmhISEkJ6ejr+/v8NRGmNMwbAEwPiUAwcO8N///peWLVuyfft2+vfvT/ny5enWrRv9+vWjS5culC9fHsBO/saYYs0SAFOspaam8umnn7Js2TLef/99AO6++25Wr15NdHQ027Zto0GDBgQGBjocqTHGFC5LAEyxc+jQIapXr05AQACTJ0/mpZdeAqBChQqMGzeOvn37epZt3LixU2EaY4yjLAEwXu/EiRMsWbKE2NhY1q1bx7lz51ixYgW/+93v6N+/PzVq1OCBBx6gYsWKTodqjDFFhuMJgIhMA+4GUoH9wFBV/cVdNwH4Pa7pgJ9Q1U/d5c35bTrgdcBo9aUhDX1cWloae/fupWHDhqgqlSpV8tR17NiRe+65h8jISACaNWtGs2bNHIrUGGOKLscTAOAzYIKqpovIy8AEYJyIhAIDgMZADeBzEWmgqhnAHGA48A2uBKAr8LEj0ZsCl56ezgcffMDXX39NTEwMx44dIzg4mG3btnHDDTfw+uuv4+/vb4PxGGNMPjieAKjqP7K8/Qbo537dC1iqqr8C8SKyD2ghIglAeVX9N4CILAR6YwlAsXH69Gk+/vhjgoKCiI6O5uuvv2bw4MGAq2f+H/7wByIiIrjQ6PPYY485Ga4xxnglxxOASzwMfOh+XRNXQnDBYXdZmvv1peXGi3355ZcsXLiQH374gU2bNgHQuXNnPvvsMyIiIpg/fz69evWiatWqDkdqjDHFQ6EkACLyOXBzDlWTVHWVe5lJQDqw+MJqOSyvlynPbd/Dcd0uoHbt2vmI2hSEs2fP8uWXX7Jy5UqSk5N59NFHad++PWvXruXtt9+mXbt2dOvWjf79+3P33XcDrt77jzzyiMORG2NM8VIoCYCqdr5cvYgMBnoCd2TpzHcYqJVlsSDgiLs8KIfy3PY9H5gPrumA8x28uWopKSns27eP8PBwEhMTqVq1Kn5+fmRmZgIQFBTEwIEDAZgwYQJTpkyhRIkSToZsjDE+w/FbACLSFRgHdFDVlCxVq4EPRGQGrk6A9YHvVDVDRJJFpBXwLfAQ8Fphx22y+/LLL3nnnXcQEbZs2cKWLVto3Lgx27dvp1SpUjzyyCNUq1aNiIgIIiMjueWWWzzrVqhQwbnAjTHGBzmeAACvAyWBz0QE4BtVHaGqO0RkGbAT162Bx9xPAACM5LfHAD/GOgAWmuTkZLZt28aOHTvYsmULS5Ys4Y033mDAgAHExcXx3nvvUblyZSpWrMjAgQM9zfjlypVj/vz5DkdvjDHmAscTAFW95TJ1U4GpOZTHAWEFGZevS0xMZOfOnSxatIiDBw8ya9YsKleuzK233uqZNOeCMmXKAPDwww9zzz33UK9ePdzJnDHGmCLK8QTAOCMzM5OjR49y6NAhfvzxR+rVq0ebNm1YuXIlI0eO5NixY55lb7zxRipUqEDlypX5+9//zoYNG+jZsydhYWEEBPz2FapQoYI15RtjjJewBKCYysjIID4+nn/961+oKiVKlPB0uHv66aeZN28eKSm/dbkYPXo0bdq04fbbb6dnz55UqFCBVq1aUa9ePcLCwjyd89q0aUObNm0c+UzGGGOuH/GlEXSjoqI0Li7O6TCumqoiIqSnp7N+/Xr27dvHgQMHCAgIYO3atcTExFC3bl3GjBnD9OnTL1o3PDycrVu38vPPPzN8+HBuueUW6tevT7Vq1WjYsCEhISGULVvWoU9mjDGmIIjIZlWNyqnOWgAKWGZmJklJSZw/f560tDSqVKlCuXLlSEhI4Mcff+T48eMkJyeTlpbGmTNnGDt2LCLCK6+8wieffEJiYiInT57k0KFD1KlTh/3795OamkrXrl2z7evYsWPUrVuXoUOHcuTIESpXrkzDhg1p27YtwcHBANSsWZO1a9cW8lEwxhhT1FgCcJXS09OpXbs2mZmZZGRkkJmZSWZmJo8//jgvvvgi8fHx1K9fn4yMjIvW++CDDxg4cCCxsbE8/PDD2bY7btw4APbt20dSUhI33HADjRo14syZM3Tq1AmA0qVLs3DhQqpXr05oaChVqlS5aD770NBQFi9enG3bxhhjzAWWAFwlPz8/evbsiZ+fH35+fvj7++Pn50dUlKulpWLFiowfP57AwECOHz9OeHg4JUqUoFWrVoBrlrpPP/0UEaFatWpUqlTJ05seYPbs2bnuW0R48MEHC/YDGmOMKdasD4AxxhhTTF2uD4BfYQdjjDHGGOdZAmCMMcb4IEsAjDHGGB/kU30ARCQROHAdN1kZOH4dt+eL7BheOzuG14cdx2tnx/DaXe9jWEdVq+RU4VMJwPUmInG5da4weWPH8NrZMbw+7DheOzuG164wj6HdAjDGGGN8kCUAxhhjjA+yBODa2AT3186O4bWzY3h92HG8dnYMr12hHUPrA2CMMcb4IGsBMMYYY3yQJQBXQUS6isgeEdknIuOdjscbiUgtEYkVkV0iskNERjsdk7cSEX8R+V5E1jgdizcSkQoiskJEdru/j62djsnbiMhT7t/j7SKyRERKOR2TNxCRd0TkmIhsz1J2k4h8JiI/uv+vWFD7twQgn0TEH5gNdANCgYEiEupsVF4pHfijqjYCWgGP2XG8aqOBXU4H4cVeBT5R1YZABHYs80VEagJPAFGqGgb4AwOcjcprLAAundt9PLBeVesD693vC4QlAPnXAtinqj+paiqwFOjlcExeR1WPqup/3K+Tcf3RrelsVN5HRIKAHsBbTsfijUSkPNAeeBtAVVNV9RdHg/JOAUBpEQkAygBHHI7HK6jql8CJS4p7Ae+5X78H9C6o/VsCkH81gUNZ3h/GTlzXRESCgabAtw6H4o1mAmOBTIfj8FYhQCLwrvs2ylsiUtbpoLyJqv4MTAcOAkeBU6r6D2ej8mrVVPUouC6UgKoFtSNLAPJPciizRymukojcAPwNeFJVTzsdjzcRkZ7AMVXd7HQsXiwAaAbMUdWmwFkKsMm1OHLfo+4F1AVqAGVF5AFnozJ5YQlA/h0GamV5H4Q1d10VEQnEdfJfrKornY7HC0UD94hIAq5bUbeLyPvOhuR1DgOHVfVC69MKXAmBybvOQLyqJqpqGrASaONwTN7sfyJSHcD9/7GC2pElAPm3CagvInVFpASuzi6rHY7J64iI4LrvuktVZzgdjzdS1QmqGqSqwbi+hxtU1a688kFV/wscEpFb3UV3ADsdDMkbHQRaiUgZ9+/1HVhHymuxGhjsfj0YWFVQOwooqA0XV6qaLiKjgE9x9XZ9R1V3OByWN4oGHgS2icgWd9lEVV3nXEjGRz0OLHYn9D8BQx2Ox6uo6rcisgL4D66ne77HRgTMExFZAnQEKovIYWAy8BKwTER+jyu56l9g+7eRAI0xxhjfY7cAjDHGGB9kCYAxxhjjgywBMMYYY3yQJQDGGGOMD7IEwBhjjPFBlgAYY4wxPsgSAGOMMcYHWQJgjLluRKS0iPxTRPxFJDjrPOf53E4JEfnSPbucMaYAWAJgjLmeHgZWqmrGtWzEPdX2euC+6xKVMSYbSwCMMXkiIrEi0sX9+k8iMiuHxQaRw9jlIhLinm73NnfLwG731LvbRWSxiHQWka9F5EcRaeFeLca9PWNMAbAEwBiTV5OBSSIyCGgKPJW10j2WfoiqJlxSfiuuWR+Hquomd/EtwKtAONAQuB9oCzwDTHQvsx24rUA+iTHGJgMyxuSNqn7pnu3taaBjDs38lYFfLimrgqtF4HeXTJoVr6rbAERkB7BeVVVEtgHB7v1liEiqiJRT1eTr/4mM8W3WAmCMyRMRaQJUB37N5YR8Dih1Sdkp4BCu2R+z+jXL68ws7zO5+MKkJHD+amM2xuTOEgBjzBWJSHVgMdALOCsid126jKqeBPxFJGsSkAr0Bh4Skfvzuc9KQKKqpl114MaYXFkCYIy5LBEpA6wE/qiqu4ApwPO5LP4PXPfyPVT1LNATeEpEeuVj152AdfkO2BiTJ6KqTsdgjCkmRKQp8LSqPngdtrUSmKCqe649MmPMpawFwBhz3ajq90CsiPhfy3bcTxTE2MnfmIJjLQDGGGOMD7IWAGOMMcYHWQJgjDHG+CBLAIwxxhgfZAmAMcYY44MsATDGGGN8kCUAxhhjjA+yBMAYY4zxQf8f3OaNcuak6SkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dpar = 1e-4\n",
    "Kp = np.hstack([(K1+K1*dpar)*np.ones([1,round(2*L/3)]), K2*np.ones([1,round(L/3)])])\n",
    "M1 = M0\n",
    "flopy.modflow.ModflowLpf(M1, hk=Kp, vka=-999., ss=-999., sy=-999., ipakcb=53)\n",
    "M1.write_input()\n",
    "success, buff = M1.run_model(silent=True)\n",
    "H1 = flopy.utils.binaryfile.HeadFile('model.hds').get_data()[0,0,:]\n",
    "dhdK = (H1-H0)/(K1*dpar)\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "plt.subplots(figsize=[8,2])\n",
    "plt.plot(X/1000., dhdK, 'k--', mfc='none', label='perturbation sensitivity, numerical')\n",
    "plt.plot(ocol/1000., dhdK[ocol], 'ro')\n",
    "plt.xlabel('$x$ (km)')\n",
    "plt.ylabel(r'$\\partial h / \\partial K$')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src onerror=\"var cell = this.closest('.code_cell');var editor = cell.querySelector('.input_area');editor.style.background='rgba(0, 200, 0, 0.2)';this.parentNode.removeChild(this)\">"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-97.274780 \n"
     ]
    }
   ],
   "source": [
    "set_background('rgba(0, 200, 0, 0.2)')\n",
    "#print('%.6f (%.6f)'% (dhdK[ocol], benchmark))\n",
    "print('%.6f '% dhdK[ocol])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Adjoint sensitivity\n",
    "\n",
    "<br>\n",
    "    \n",
    "$\\begin{align}\n",
    "\\quad\\quad \\frac{\\partial h(x')}{\\partial K}=\\int\\limits_X\\psi^*_1(x)\\ b\\ \\frac{d^2 h(x)}{d x^2}\\; dx\n",
    "\\cr\\end{align}$\n",
    "\n",
    "Governing equation:\n",
    "\n",
    "$\\begin{align}\n",
    "\\quad\\quad K\\ b\\ \\frac{d\\psi^*_1}{dx}+\\frac{1}{2\\ K\\ b}\\delta(x-x')=0\\cr\n",
    "\\cr\\end{align}$\n",
    "\n",
    "Boundary conditions: \n",
    "\n",
    "$\\begin{align}\n",
    "\\quad\\quad -K\\ b\\ \\frac{d\\psi^*_1(x)}{dx}=0\\ ,\\quad\\quad\\; x=0=\\Gamma_2\\cr\n",
    "\\end{align}$\n",
    "\n",
    "$\\begin{align}\n",
    "\\quad\\quad \\psi^*_1(x)=0\\ ,\\ \\quad\\quad\\quad\\quad\\quad x=L=\\Gamma_1\\cr\n",
    "\\cr\\end{align}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3a. Analytical\n",
    "\n",
    "$\\begin{align}\n",
    "\\quad\\quad\\rm{Not\\ available}\n",
    "\\cr\\end{align}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3b. Semi-analytical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfsAAACdCAYAAACzUdEhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjf0lEQVR4nO3deXwV5dn/8c+VhD0poKwCNaFFrRIIiwHBCgqKiIjSuvBgQdT6QFvqDx5RlKIU10rrxkulFln0sVVBVKAKohJwQRQpLuyLAQJKImtAQxau3x85nCeHBAiQ5OScfN+vV16Zue85M9eZk5Nr7pl75jZ3R0RERKJXTLgDEBERkfKlZC8iIhLllOxFRESinJK9iIhIlFOyFxERiXJK9iIiIlEuLtwBlJcGDRp4YmJiuMMQERGpMJ9//vn37t7wyPKoTfaJiYksW7Ys3GGIiIhUGDPbXFK5TuOLiIhEOSV7ERGRKBe1p/HL2sSJE5kzZ05IWfXq1Zk7dy4Ajz76KO+9915Ifd26dXn11VcBGD9+PB9//HFIfZMmTZg2bRoAd999NytWrAipT0pK4plnngFgxIgRrF27NqT+3HPP5a9//SsAw4YNY8uWLSH1HTp0YPz48QAMGTKErKyskPoLL7yQ0aNHA3DDDTdw4MCBkPqePXty++23A3D11Vdz6NChkPq+ffvy29/+loMHD3LttddypOuuu44bb7yRvXv3MmjQoGL1gwYN4le/+hXfffcdQ4cOLVZ/2223ccUVV/DNN98wcuTIYvV//OMfufjii1m9ejVjxowpVn/XXXfRqVMnli9fzoMPPlis/r777qNNmzZ8/PHHPPbYY8XqH374YVq1asX777/Ps88+W6z+scceo0WLFrz11lvBz7GoZ555hgYNGjBr1ixeeeWVYvXPP/888fHx/POf/2T27NnF6l966SViY2OZMmUKCxYsCKmrXr0606dPD27nww8/DKn/yU9+wqRJkwB4/PHHi13SatSoEY8//njwfa5cuTKk/qc//SkPPfQQAOPGjWPjxo0h9WeddRZjx44FYPTo0Wzfvj1YV6NGDQYOHEj37t2LvScRCQ8l+1LKyclh//79IWXVq1cPqd+3b19IfWxsbHD6wIED7NmzJ6S+du3awens7Gx27doVUn/aaacFp/fs2cP3338fUr93797g9M6dO8nMzAypL7q9zMxMduzYcdT67777juzs7KPWb9u2jYKCghLr3Z2MjAyOdDi+goKCYgciQHB/5efnk56eXqz+cDy5ubls2rSpWP3hz+PHH39kw4YNR60/cOAA69atK1Z/+OBm3759rFmzplj9jz/+CBS+z1WrVhWrP3jwIFC477/++uti9bm5uQDs2LGDL7/8slj94f25ffv2Ygd6ULhfAbZs2cLy5ctD6mrUqBGcTk9PL5bMi/7trF+/nk8//TSkvnnz5sHp1atX88knn4TU7969Ozj95Zdf8sUXX4TUH37vAMuXLw/Z/7t27WLfvn1K9iKViEXrQDgdO3Z0ddATqXg5OTns3buXxo0bs3LlSm6++WbGjBlD3759MbNwhycS1czsc3fveGS5rtmLSJmqWbMmjRs3BgrPKH3//ff069ePlJQUZsyYUexykIiUPyV7ESk3F198MWvXrmX69Onk5ORw3XXXkZqaqoQvUsGU7EWkXMXFxTFo0CBWrVrFyy+/zMCBA4mJicHdef3118nLywt3iCJRTx30RKRCxMbGcv311wfnP/74Y/r378+ZZ57J6NGjGTJkSEjHQxEpO2rZi0hYdOnShblz59KkSROGDRtGy5YtefLJJ8nJyQl3aCJRR8leRMLCzOjTpw9LlixhwYIF/PznP+fhhx8O3nIYrXcKiYSDkr2IhJWZ0bNnTxYtWsSKFSuoVasWeXl5dOrUiQceeCDkeRIicnKU7EWk0mjSpAlQ+FCfxo0bM3bsWM4880zuvfdedu7cGeboRCKXkr2IVDqNGjVizpw5LF++nB49enD//feTmJhY7LG+IlI6SvYiUmm1a9eO1157ja+//pqhQ4fyi1/8AoB58+axbdu2MEcnEjmU7EWk0jvvvPOYMGECMTEx5ObmMmjQIFq2bMmwYcPYvLnE4btFpAglexGJKNWrV2fp0qUMGTKE559/np///OfccsstJQ6mJCKFlOxFJOIkJSUxadIkNm3axLBhw/jnP/8ZHNVRj+IVKU7JXkQiVvPmzXnqqafYtm0bnTp1AmDYsGFcd911xYblFanKwpbszSzWzP5jZnMD86eZ2QIzWx/4Xb/Isneb2QYzW2tmvcIVs4hUTqeddlpwulmzZsyfP5+UlBT69evHZ599FsbIRCqHcLbsbwdWF5kfDbzn7q2A9wLzmNm5wA3AecDlwDNmFlvBsYpIhLj33nvZvHkz48eP54MPPiA1NZUnn3wy3GGJhFVYkr2ZNQf6AJOLFPcDpgempwNXFyl/2d0Puvs3wAYgtYJCFZEIVK9ePcaOHcvmzZt55JFHuOqqqwD44osveP/99/UoXqlywtWyfwK4Eyjak6axu38LEPjdKFDeDNhaZLmMQJmIyDElJCRw1113kZSUBMCECRPo0aMHF154IfPmzVPSlyqjwpO9mV0JZLr756V9SQllJX5Dzew2M1tmZsuysrJOOkYRiU6TJ0/m6aefZuvWrfTu3ZvU1FTmzZsX7rBEyl04WvZdgavMLB14GbjEzP4X2GFmTQECvzMDy2cALYq8vjmwvaQVu/tz7t7R3Ts2bNiwvOIXkQhVs2ZNfve737FhwwYmT57Mrl27gh34Dh06REFBQZgjFCkfFZ7s3f1ud2/u7okUdrx7391vBGYDgwOLDQbeDEzPBm4wsxpmlgS0Aj6t4LBFJIpUr16dW265hbVr13LHHXcAMGvWLFq3bs2LL75Ifn5+mCMUKVuV6T77R4BLzWw9cGlgHndfCbwKrALmAb93dx1+i8gpi4uLo1atWgDUrVuXatWqMWjQIM455xyef/55cnNzwxyhSNmwaO2g0rFjR1+2bFm4wxCRCHLo0CHmzJnD/fffz+eff84ll1zCe++9F+6wRErNzD53945HlseFIxgRkcooJiaGfv36cdVVVzFv3rzgo3d/+OEHpkyZws0330zt2rXDHKXIiatMp/FFRCoFM6N379706dMHgDlz5jB8+HASExP5y1/+QnZ2dpgjFDkxSvYiIsdx/fXXs3jxYtq1a8fo0aNJTEzk/vvvV0c+iRgnlezNrI4eWSsiVckvf/lL5s+fz9KlS+natSvz588nNrbw3+DBgwfDHJ3IsZUq2ZtZjJn9l5n928wygTXAt2a20swmmFmr8g1TRKRySE1NZfbs2bz77ruYGd999x0tWrTgzjvvDA6zK1LZlLZlvxD4GXA30MTdW7h7I+CXwCfAI2Z2YznFKCJS6dSsWROA/Px8Lr30Uv72t7+RmJjI7bffzrZt28IcnUioUt16Z2bV3D3vVJepSLr1TkQq0vr163n44Yd58cUXiYuLY8uWLehJnlLRjnbrXala9qVJ4pUp0YuIVLRWrVoxZcoU1q9fz8SJE5XopVI5brI3s0vN7B9mlhKYv63coxIRiVCJiYnceuutAPTq1YsBAwaEOSKR0j1U53fAEOBPZnYakFKuEYmIRIndu3cTE6M7nCX8SvNXmOXue9z9DuAy4PxyjklEJCqYlTRCt0jFK02y//fhCXcfDbxQfuGIiESXaB1/RCLLcZO9u795xPzE8gtHRCR6qGUvlcUpD4RjZh+5e9eyCEZEJJpcfvnlGjhHKoWyGPXujDJYh4hI1Bk3bly4QxABSpnszWwi8FXg52t3Lzrkky5IiYiIVGKlbdl/BbQBBgKtzWwf/5f8E8opNhGRiHbxxRdTp04d5s6dG+5QpIorVbJ39+eKzptZcwqTfzIwvxziEhGJeAcPHiQuriyuloqcmtKexr8A+MQD95C4ewaQAbxVjrGJiEQ09caXyqK0j3YaDCw3s5fN7CYza1KeQYmIRAvdZy+VQWlP4w8FMLNzgN7ANDOrS+HQt/OAj9y9oNyiFBGJQGrZS2VR2tP4ce6e7+5rgDXA42ZWC7gYuBZ4DCg2pJ6ISFV2zTXX6Jq9VAqlHc9+ubu3r4B4yozGsxcRkarmlMazB3QuSkTkBBUUFFBQoCucEn6lPb/U0MxGHq3S3R8ro3hERKLGJZdcQkxMDAsXLgx3KFLFlTbZxwLxqIUvIiIScUqb7L919/HlGomISBTSrXdSGeiavYhIOdGtd1JZlDbZX3K8BUx/1SIiIpVSaU/jv25mrwFvuvuWw4VmVh24kMIn7C0EppV5hCIiEWrAgAE6jS+VQmmT/eXAzcC/zCwJ2APUpLDj3jvA4+6+orQbNbMWwAtAE+AQ8Jy7P2lmpwGvAIlAOnCdu+8OvOZu4BagAPiju2sAHhGp1P77v/873CGIAKV/XG4O8AzwjJlVAxoAP7r7npPcbj7wP+6+3MwSgM/NbAFwE/Ceuz9iZqOB0cBdZnYucANwHnAG8K6ZnaVH9IpIZbZ//37cnYQEjQQu4VXaa/ZB7p7n7t+eQqIn8PrlgelsYDXQDOgHTA8sNh24OjDdD3jZ3Q+6+zfABiD1ZLcvIlIR+vXrxxVXXBHuMEROPNmXNTNLBNoBS4HG7v4tFB4QAI0CizUDthZ5WUagTESkUtM1e6kMwprszSweeA34f+6+71iLllBW7BtkZreZ2TIzW5aVlVVWYYqInBTdpCSVRdiSfeDa/2vAS+4+K1C8w8yaBuqbApmB8gygRZGXNwe2H7lOd3/O3Tu6e8eGDRuWX/AiIqWklr1UBsdN9mZWJ/A7vqw2Grgn/3lg9RHP1Z9N4W18BH6/WaT8BjOrEbgboBXwaVnFIyJSHtSyl8qiNL3x65vZEAo7xc0ro+12BX4DfGVmKwJl9wCPAK+a2S3AFuBaAHdfaWavAqso7Mn/e/XEF5HKbvDgweTl5YU7DJHjj2dvZoOB4cAUYKa7Zx5j2bvc/S9lG+LJ0Xj2IiJS1RxtPPvStOw/pfCBOi2OTPSB1nZwFkgBKkWyFxEJt6ysLA4dOkTjxo3DHYpUccdN9u6+OjD5ZQnV+9z91sMzZvZsWQUmIhLpbrzxRvbt28eSJUvCHYpUcaV9XC4AZvYvYAWFiX8Z8MARi4wpm7BERKKDeuNLZXCit95NAn4AfgUsAB4PPO4WAHffVYaxiYhENPXGl8riRJP9T4FFwFB3TwFmAn8u66BERKKFWvZSGZzQaXzgLApvhzvXzLIpPJ3fzcz+DXzp7npsnYhIgFr2UlmcaLJ/6nBCN7P6QEegNzCQwnvkzy/b8EREItdtt91GTk5OuMMQOeFk/04gya8G1gLnAP9y99vLPDIRkQh3zTXXhDsEEeAEk727tzOzWAqT/NnAu5TdU/VERKLKli1byM/Pp2XLluEORaq4E23ZE3hM7crAj4iIHMWwYcPIzMzks88+C3coUsWFfTx7EZFopt74Uhko2YuIlBP1xpfKQsleRKQcqWUvlYGSvYhIOVHLXiqLE+6gJyIipTN8+HB++OGHcIchomQvIlJeLrvssnCHIALoNL6ISLlZt24dX331VbjDEKlaLfu8vDwyMjL0+Eqp8mrWrEnz5s2pVq1auEOJaqNGjWLLli385z//CXcoUsVVqWSfkZFBQkICiYmJ6jgjVZa7s3PnTjIyMkhKSgp3OFFPvfGlMqhSp/FzcnI4/fTTleilSjMzTj/9dJ3hqgD6XyOVRZVK9qAvnwjoe1CR1LKXyqDKJftIMG3aNP7whz8AMGnSJF544YVjLn/vvffy7rvvHnOZtLQ0Pv7441OOLT09ndatWx93uYceeihkvkuXLqe87YpQmv19qop+vsdaZvv27cH5W2+9lVWrVp3wttLS0rjyyitP+HVSNnRQJZVFlbpmH4mGDh163GXGjx9/3GXS0tKIj4+vsKT70EMPcc899wTny+JAoyKUZn9XhGnTptG6dWvOOOMMACZPnhzmiORk3HHHHezfvz/cYYioZV/Rrr76ajp06MB5553Hc889FyyfOnUqZ511Ft26deOjjz4Klo8bN46//vWvAKxYsYLOnTvTpk0brrnmGnbv3g3ATTfdxMyZMwFITEzkvvvuo3379iQnJ7NmzRrS09OZNGkSjz/+OCkpKXzwwQchMS1atIiUlBRSUlJo164d2dnZuDujRo2idevWJCcn88orrxR7L0e2UK+88krS0tIYPXo0P/74IykpKQwcOBCA+Ph4gKOuNy0tje7du/PrX/+ac845h4EDBx739OeMGTNo3bo1bdu25aKLLgKgoKCAUaNGcf7559OmTRv+/ve/B9ffrVs3rrvuOs466yxGjx7NSy+9RGpqKsnJyWzcuLHY/i7tZxcfH8+YMWNo27YtnTt3ZseOHQDMmTOHTp060a5dO3r27BksPyw7O5ukpCTy8vIA2LdvH4mJicyYMYNly5YxcOBAUlJS+PHHH+nevTvLli0DYN68ebRv3562bdvSo0cPAD799FO6dOlCu3bt6NKlC2vXrj3mvpOK0bVrV3r16hXuMEQK//lG40+HDh38SKtWrQqZ79atW7Gfp59+2t3dDxw4UGL91KlT3d09KyurWF1p7Ny5093df/jhBz/vvPP8+++/9+3bt3uLFi08MzPTDx486F26dPHf//737u5+3333+YQJE9zdPTk52dPS0tzdfezYsX777be7u/vgwYN9xowZ7u5+5pln+lNPPeXu7k8//bTfcsstxdZzpCuvvNI//PBDd3fPzs72vLw8nzlzpvfs2dPz8/P9u+++8xYtWvj27dv9m2++8fPOO8/d3adOnRqM0929T58+vnDhQnd3r1OnTsg2Ds8fbb0LFy70n/zkJ75161YvKCjwzp07+wcffHDMfdm6dWvPyMhwd/fdu3e7u/vf//53v//++93dPScnxzt06OCbNm3yhQsXet26dX379u2ek5PjZ5xxht97773u7v7EE08E9+Wx9lNJn527O+CzZ892d/dRo0YFt79r1y4/dOiQu7v/4x//8JEjRxbbbzfddJO//vrrwdgPL9OtWzf/7LPPgts+PJ+ZmenNmzf3TZs2hcS0d+9ez8vLc3f3BQsWeP/+/d3dfeHChd6nT58S38+R3wcpe1988YUvWbIk3GFIFQIs8xJyolr2Feypp54KtgC3bt3K+vXrWbp0Kd27d6dhw4ZUr16d66+/vtjr9u7dy549e+jWrRsAgwcPZvHixSVuo3///gB06NCB9PT048bUtWtXRo4cyVNPPcWePXuIi4vjww8/ZMCAAcTGxtK4cWO6detWJmNyH2u9qampNG/enJiYGFJSUo4be9euXbnpppv4xz/+QUFBAQDvvPMOL7zwAikpKXTq1ImdO3eyfv16AM4//3yaNm1KjRo1+NnPfhZ8ullycnKp9lNJnx1A9erVg9fFi+7zjIwMevXqRXJyMhMmTGDlypXF1nnrrbcydepUoPDszpAhQ44ZwyeffMJFF10UvGXutNNOAwr/Pq699lpat27NiBEjStyWVLw///nP/Pa3vw13GCJV+5p9WlraUetq1659zPoGDRocs/5o23v33XdZsmQJtWvXpnv37sHbn8qyI0+NGjUAiI2NJT8//7jLjx49mj59+vDWW2/RuXNn3n333VL1II6Li+PQoUPB+dLcynWs9R6OG0oX+6RJk1i6dCn//ve/SUlJYcWKFbg7EydOLHbqNC0tLWT9MTExwfmYmJhi29q6dSt9+/YFCq/jn3POOUf97KpVqxb8/IrGPXz4cEaOHMlVV11FWloa48aNK/YeunbtSnp6OosWLaKgoOC4nR/dvcS/lbFjx3LxxRfz+uuvk56eTvfu3Y+5Hqk4pfkuiZQ3tewr0N69e6lfvz61a9dmzZo1fPLJJwB06tSJtLQ0du7cSV5eHjNmzCj22rp161K/fv3g9fYXX3wx2MovjYSEBLKzs0us27hxI8nJydx111107NiRNWvWcNFFF/HKK69QUFBAVlYWixcvJjU1NeR1iYmJrFixgkOHDrF161Y+/fTTYF21atWC16KLKs16j3T33Xfz+uuvlxh3p06dGD9+PA0aNGDr1q306tWLZ599NrjtdevWceDAgePunyO1aNGCFStWsGLFCoYOHXrUz+5Y9u7dS7NmzQCYPn36UZcbNGgQAwYMCGnVH+3zuuCCC1i0aBHffPMNALt27Sq2rWnTppX6fUr5MjO+//57Jk2aFOyot3nzZpYuXcqqVavYsmULu3fvLvG7IlKWqnTLvqJdfvnlTJo0iTZt2nD22WfTuXNnAJo2bcq4ceO44IILaNq0Ke3btw+elob/a/VPnz6doUOH8sMPP9CyZcvg6d/S6Nu3L7/+9a958803mThxIr/85S+DdU888QQLFy4kNjaWc889l969e1O9enWWLFlC27ZtMTMeffRRmjRpEnK6u2vXriQlJZGcnEzr1q1p3759sO62226jTZs2tG/fnpdeeilYfs0115S43jVr1hw19q+++oqrrrqqWPmoUaNYv3497k6PHj1o27Ytbdq0IT09nfbt2+PuNGzYkDfeeKPU++lojvbZHcu4ceO49tpradasGZ07dw4m6CMNHDiQP/3pTwwYMCBYdtNNNzF06FBq1arFkiVLguUNGzbkueeeo3///hw6dIhGjRqxYMEC7rzzTgYPHsxjjz3GJZdccsrvV8pGq1ateO211xg2bBj9+/cnPj6eyZMn88ADDxRbNjs7m/j4eB588EFefvllEhISiI+PJz4+noSEBKZNm4aZ8fbbb7Nu3bpgeXx8PPXq1QveabNv3z5iY2OpXbu2bv2TIIvWU0wdO3b0w72XD1u9ejW/+MUvwhTRyRk+fDjt27c/7rXcaNarVy/mz58f7jDKzcyZM3nzzTd58cUXK3S7kfh9iET79+9n//79NGrUiJiYGDZu3MjatWvZv38/2dnZwd/33HMPMTExTJ06lblz54bUFRQUBJ+zcOONN4YcQEPhZcWsrCyg8ID6jTfewMyCBwtnn302CxcuBAov+axfvz54oHD4EeI333wzAIsXLyY3NzfkYKNevXrUrVu3AveanCwz+9zdOxYrV7KvvMaOHcv8+fN5++23Of3008MdjpSD4cOH8/bbb/PWW29x1llnVei2I+37IIXy8vKKHSjk5+fTtWtXAObOncvq1atDlqlbty4TJkwA4Oabb+ajjz4KqT///PNZunQpACkpKXzxxRch2+zevXvwYKFTp07s2LEj5MxC9+7dGTNmDAAPPvgg7h5yVqJVq1bBM38bN26kdu3aJCQkULt2bWJidDW5LEV8sjezy4EngVhgsrs/cqzloyHZi5QnfR8ECjsQ5ubmBjusrlq1il27doUcTDRs2DDYYXXMmDFs27YtpP7CCy/k0UcfBaBRo0bBswyHDRkyhClTpuDuVKtWLXiZ0syoU6cOw4cP56GHHiI3N5dLL7005EAhPj6eK664gssuu4ycnBxeffXVkLMS8fHxNGvWjPr16x+1A2tVcrRkHxHX7M0sFngauBTIAD4zs9nufuLPDxURkSAzC7lT5dxzzz3m8g8++OAx6zMzM8nLy+PAgQPBA4I6deoAhQcW06dPDzlQOHxmASA3N5eYmBgyMzPZuHFjcJlGjRpx2WWXkZWVxeDBg4tt829/+xsjR45k3bp1JCcnF+vvMGbMGPr27Ut6ejqPPvpoyIFCQkICPXr0ICkpib1797Jx48aQ19apUycqzj5ERLIHUoEN7r4JwMxeBvoBJ5zsdeQnotvBpHxVq1aNevXqUa9evZDymJiY4FM1SxIfHx+8XFCSpk2bsmHDhpADhf3795OcnAwU3rV0xx13FDuYOHwwk5mZycyZM8nOzg65VXjWrFkkJSWxZMkSevfuXWy7CxYsoGfPnrzzzjvcddddxc4sjBkzhpYtW/L111+zaNGiYgcTKSkp1KpVi4MHDwKFz+ao6DwUKcm+GbC1yHwG0OlEV1KzZk127typYW6lSvPAePY1a9YMdygiJyQuLo6f/exnR61v0qRJsUG4ikpNTSUzMxOA/Pz84MFC/fr1AWjfvj1vvvlmsT4Rh7dZs2ZNWrRoQXZ2NllZWWzatIn9+/czYsQIoLBzY0mDXK1bt45WrVoxceJERo0aRVxcHAkJCcydO7fCxiuJiGv2ZnYt0Mvdbw3M/wZIdffhRyx3G3AbwE9/+tMOmzdvDllPXl4eGRkZGsdbqryaNWvSvHlzqlWrFu5QRKJGbm4ue/bsCR4kHD5g6NatG7Vq1WLp0qW8//77wfIRI0YEn4ZZViK6g56ZXQCMc/degfm7Adz94aO9pqQOeiIiItHsaMk+UnodfAa0MrMkM6sO3ADMDnNMIiIiESEirtm7e76Z/QGYT+Gtd1PcXSN9iIiIlEJEJHsAd38LeCvccYiIiESaiLhmfzLMLAvYfNwFS68B8H0Zrq+q0n48ddqHp0778NRpH5668tiHZ7p7wyMLozbZlzUzW1ZSpwc5MdqPp0778NRpH5467cNTV5H7MFI66ImIiMhJUrIXERGJckr2pfdcuAOIEtqPp0778NRpH5467cNTV2H7UNfsRUREopxa9iIiIlFOyb4UzOxyM1trZhvMbHS444k0ZtbCzBaa2WozW2lmt4c7pkhlZrFm9h8zmxvuWCKRmdUzs5lmtibw93hBuGOKRGY2IvBd/trM/mVmGlXpOMxsipllmtnXRcpOM7MFZrY+8Lt+eW1fyf44zCwWeBroDZwLDDCzYw/4LEfKB/7H3X8BdAZ+r3140m4HVoc7iAj2JDDP3c8B2qJ9ecLMrBnwR6Cju7em8KmmN4Q3qogwDbj8iLLRwHvu3gp4LzBfLpTsjy8V2ODum9w9F3gZ6BfmmCKKu3/r7ssD09kU/oNtFt6oIo+ZNQf6AJPDHUskMrOfABcBzwO4e6677wlrUJErDqhlZnFAbWB7mOOp9Nx9MbDriOJ+wPTA9HTg6vLavpL98TUDthaZz0CJ6qSZWSLQDlga5lAi0RPAncChMMcRqVoCWcDUwKWQyWZWJ9xBRRp33wb8FdgCfAvsdfd3whtVxGrs7t9CYaMIaFReG1KyPz4roUy3MJwEM4sHXgP+n7vvC3c8kcTMrgQy3f3zcMcSweKA9sCz7t4OOEA5njaNVoHryv2AJOAMoI6Z3RjeqOR4lOyPLwNoUWS+OTpldcLMrBqFif4ld58V7ngiUFfgKjNLp/BS0iVm9r/hDSniZAAZ7n74rNJMCpO/nJiewDfunuXuecAsoEuYY4pUO8ysKUDgd2Z5bUjJ/vg+A1qZWZKZVaewI8rsMMcUUczMKLxOutrdHwt3PJHI3e929+bunkjh3+D77q7W1Alw9++ArWZ2dqCoB7AqjCFFqi1AZzOrHfhu90AdHU/WbGBwYHow8GZ5bShihrgNF3fPN7M/APMp7HU6xd1XhjmsSNMV+A3wlZmtCJTdExi2WKQiDQdeChy4bwKGhDmeiOPuS81sJrCcwjtt/oOepndcZvYvoDvQwMwygPuAR4BXzewWCg+iri237esJeiIiItFNp/FFRESinJK9iIhIlFOyFxERiXJK9iIiIlFOyV5ERCTKKdmLiIhEOSV7ERGRKKdkLyInxcxqmdkiM4s1s8Si43Sf4Hqqm9niwAhqIlIOlOxF5GTdDMxy94JTWUlg6Oj3gOvLJCoRKUbJXkSKMbOFZnZpYPoBM3uqhMUGUsKzvM2sZWAI2fMDLf41geFkvzazl8ysp5l9ZGbrzSw18LI3AusTkXKgZC8iJbkPGGNmA4F2wIiilYFny7d09/Qjys+mcHTDIe7+WaD458CTQBvgHOC/gAuBO4B7Ast8DZxfLu9ERDQQjogU5+6LAyOajQS6l3CqvgGw54iyhhS29H91xGBR37j7VwBmthJ4z93dzL4CEgPbKzCzXDNLcPfssn9HIlWbWvYiUoyZJQNNgYNHSb4/AjWPKNsLbKVwlMOiDhaZPlRk/hChDY4aQM7JxiwiR6dkLyIhzKwp8BLQDzhgZr2OXMbddwOxZlY04ecCVwODzOy/TnCbpwNZ7p530oGLyFEp2YtIkJnVBmYB/+Puq4H7gXFHWfwdCq+9B7n7AeBKYISZ9TuBTV8MvHXCAYtIqWg8exE5KWbWDhjp7r8pg3XNAu5297WnHpmIHEktexE5Ke7+H2ChmcWeynoCPfvfUKIXKT9q2YuIiEQ5texFRESinJK9iIhIlFOyFxERiXJK9iIiIlFOyV5ERCTKKdmLiIhEOSV7ERGRKPf/AWEOL+ns4KDzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "M0 = timml.Model3D(kaq=K1, z=[0., -b])\n",
    "timml.inhomogeneity1d.StripInhom3D(M0, x1=2*L/3., x2=L, kaq=K2, z=[0,-b], topboundary='conf')\n",
    "timml.ImpLineDoublet1D(M0, xld=0.)\n",
    "timml.HeadLineSink1D(M0, xls=L, hls=0.)\n",
    "timml.LineSink1D(M0, xls=float(ocol), sigls=-1.)\n",
    "M0.solve(silent=True)\n",
    "A = M0.headalongline(X, 0.).flatten()\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "plt.subplots(figsize=[8,2])\n",
    "plt.plot(X/1000., A, 'k--', mfc='none', label='adjoint solution, semi-analytical')\n",
    "plt.xlabel('$x$ (km)')\n",
    "plt.ylabel(r'$\\psi^*_1$ (T/$L^2$)')\n",
    "plt.legend(loc=3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgAAAACdCAYAAADCH8BOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAArHElEQVR4nO3dfXzP9f748cdzGxtGdpCLDaMYdsFsiRQ6lchchKKkdOVUkpNT4RR1LpxydEpLLjs/uiBfqYNQ7ahEHSpDGZYI+TAXcznErp6/Pz6ffdpmm9nVZxfP++32ufl8Xq/P+/1+vuezvZ6f9+viLaqKMcYYY6oWL08HYIwxxpiyZwmAMcYYUwVZAmCMMcZUQZYAGGOMMVWQJQDGGGNMFWQJgDHGGFMF+Xg6gLJUv359DQ4O9nQYxhhjTJmIj49PVtUGedVV6ARARHoBrwHewJuq+lJB7w8ODmbjxo1lEpsxxhjjaSKyL7+6CtsFICLewBtAb6AdcJeItPNsVMYYY0zFUGETAKATsEtVf1bVVGAR0N/DMRljjDEVQkXuAggE9md77QCuzf0mERkJjARo1qxZ2URmjDHlnMPh4Pz586gqqoqfn5/7b+SePXv49ddf3XWqSq1atWjZsiUA27dvz7GtqlK3bl1atWoFQHx8PKmpqTnqGzRoQEhICABr164lIyMjR31gYCBt27ZFVYmLi8tRp6q0bNmStm3bkpaWxsqVKwFy1Ldt25Z27dpx9uxZli5dSmZmpvuRkZHBtddeS3h4OMePH2fhwoU56jMzM+nZsycREREcOHCAefPm5dg2MzOTIUOGEBERwU8//cSsWbMu2v7RRx8lLCyMzZs3M336dPd2WY9JkybRpk0bvvzyS1577bWLtn/99de56qqryvQzUJETAMmj7KIbG6jqHGAOQHR09EX1aWlp7l8EY0zh+Pn5ERQURLVq1TwdiimCP//5z7z44os5yqKjo/nuu+8AGDhwIFu2bMlRf+ONN/L5558D0K9fP3bv3p2jvl+/fixbtgyA2267jSNHjuSoHzZsGO+++y4AvXr14tdff81R/4c//IFZs2ahqvTq1euimJ966immTp3KuXPnuP322y+q/8tf/sKkSZM4ceIE99xzz0X1r7zyCuHh4Rw+fJjRo0dfVF+3bl13AjBx4kR3uYjg7e1N+/btiYiI4ODBg8yZMwcvL68cjwEDBhAWFsbRo0eJi4tzl3t7e+Pl5UVKSgoAKSkp7Nq166Lt09LSLoqptFXkBMABNM32Ogg4eNk7cTioXbs2wcHBiOSVUxhjslNVjh07hsPhoEWLFp4OxxRBaGgoDzzwAD169EBEEBHq1avnrn/ppZc4efKku05EuPLKK931M2fO5Ny5cznqGzZs6K5/7733SE1NzVHfuHFjd/0nn3xCZmZmnvUiwv/+9z/389z1/v7+bNq0yf33Oqu+UaNGADRq1IidO3de1MDWrVsXgFatWnH06NGL6n19fQFnIpSWloaXl5d739l1797d3ZjnpWfPnuzfvz/f+piYGGJiYvKtL0tSUe8GKCI+wE7gJuAA8B1wt6puy2+b6OhozT0LYMeOHbRp08Yaf2Mug6qSmJhI27ZtPR2KuQxnzpzB39/f02GYMiQi8aoanVddhR0EqKrpwOPAp8AOYHFBjX9BrPE35vLY70zFk5CQQMuWLfnoo488HYopJypsAgCgqqtUtbWqXqWqkz0dT2maP38+jz/+OACzZs3i7bffLvD9kyZNYvXq1QW+Z82aNe5LbWVp+fLlvPSSc8mGpUuXsn37dnddYeK+7bbbOHnyJCdPnmTGjBmlGmtRXHfddaV+jODgYJKTk/Ot37t3LwsXLnS/3rhxI0888USRjtWjRw9bP6OCO3ToEDExMfj4+NChQwdPh2PKiYo8BqDKeuSRRy75nr/+9a+XfM+aNWvw9/cvkwYru379+tGvXz/AmQDExMTQrp1zCYfCxL1q1SrA2cjNmDGDxx57rPSCLQJPJFW5ZSUAd999N+Ds14yOzvMqoKnkzp07R79+/Th69Chr166ladOml97IVAkV+gpAZTFgwACioqIIDQ1lzpw57vJ58+bRunVrunfvztdff+0uf+GFF3j55ZcB2LJlC507dyYiIoLbb7+dEydOADBixAiWLFkCOL8tPv/883Ts2JHw8HASExPZu3cvs2bN4tVXX6VDhw6sW7cuR0xffvklHTp0oEOHDkRGRroHvUydOpVrrrmGiIgInn/+ecDZ2LRt25aHH36Y0NBQevbs6R7hGxsbS7t27YiIiGDo0KHAb1cz/ve//7F8+XKefvppOnTowO7du91xf/zxx9x5553ueNasWUPfvn3d55OcnMz48ePZvXs3HTp04Omnn2b48OHuUcjgHHW8fPnyfH/uZ8+epU+fPrRv356wsDD+7//+D3BOYerevTtRUVHceuutJCUlAc5vwk8++STdunWjbdu2fPfddwwcOJBWrVrx3HPPufebXx/rRx99xLXXXktkZCQ333wzhw8fdv9/Zg3IatmyJbGxsZf8bGSZOHEir732mvv1s88+S2xsLOPHj2fdunV06NCBV199lTVr1rgHHp05c4b777+f8PBwIiIi+OCDDwB49NFHiY6OJjQ01P1/ayq2zMxMhg8fzsaNG1m4cCFRUVGeDsmUJ7nnWlbmR1RUlOa2ffv2i8rK2rFjx1RV9dy5cxoaGqrJycl68OBBbdq0qR45ckQvXLig1113nY4aNUpVVZ9//nmdOnWqqqqGh4frmjVrVFV14sSJOmbMGFVVve+++/T9999XVdXmzZtrbGysqqq+8cYb+uCDD160n9xiYmL0q6++UlXVlJQUTUtL008//VQffvhhzczM1IyMDO3Tp49++eWXumfPHvX29tbNmzerquodd9yh77zzjqqqNm7cWM+fP6+qqidOnFBV1Xnz5rnPJXuc2V+npaVp06ZN9cyZM6qq+sgjj7j32bx5cz169Kju2bNHQ0ND3duuWbNG+/fvr6qqJ0+e1ODgYE1LS8v3575kyRJ96KGH3K9Pnjypqamp2qVLFz1y5Iiqqi5atEjvv/9+VVXt3r27PvPMM6qqOm3aNG3cuLEePHhQz58/r4GBgZqcnKyqqrVq1crzeMePH9fMzExVVZ07d66OHTtWVZ3/D126dNHz58/r0aNH9Xe/+52mpqaqat6fjdw/g8jISFVVzcjI0JYtW2pycrJ+8cUX2qdPH/exs79+5pln3J+TrLiyHys9PV27d++u33//vfu8v/vuu4vOpzz87piCZWRk6Lhx4/SVV17xdCjGQ4CNmk+baF0AufTo0eOisjvvvJPHHnuMc+fOcdttt11UP2LECEaMGEFycjKDBw/OUbdmzZpLHjM2Npb//Oc/AOzfv5+ffvqJQ4cO0aNHDxo0cN7DYciQIezcuTPHdqdOneLkyZN0794dgPvuu4877rgjz2MMHDgQgKioKD788MNLxtS1a1fGjh3LsGHDGDhwIEFBQcTFxREXF0dkZCTg/Cb5008/0axZM1q0aOHuW4yKimLv3r0AREREMGzYMAYMGMCAAQMuedwsPj4+9OrVi48++ojBgwezcuVK/vnPfxa4Tffu3Rk1ahRHjhzhww8/ZNCgQfj45P8RDw8P56mnnmLcuHHExMRwww03kJCQQEJCArfccgsAGRkZOaYvZXVdhIeHExoa6q5r2bIl+/fvzzGVKjeHw8GQIUNISkoiNTU1xxS6Pn364Ovri6+vL1deeSWHDx8mKCgoz89G9mMEBwdTr149Nm/ezOHDh4mMjCwwBoDVq1ezaNEi9+uAgAAAFi9ezJw5c0hPTycpKYnt27cTERFR4L5M+XX+/Hn8/Pzc422Myc26ADxszZo1rF69mvXr1/P9998TGRnpXpSoJEdaZ81x9fb2Jj09/ZLvHz9+PG+++Sa//vornTt3JjExEVVlwoQJbNmyhS1btrBr1y4efPDBHPvPfYyVK1cyatQo4uPjiYqKKtSxswwZMoTFixfz+eefc80111C7du1LbjN8+HAWLFjAvHnzuP/++wt8b+vWrYmPjyc8PJwJEybw17/+FVUlNDTUfY5bt24lLi7OvU3WeWafN5z1Ove5Pfvss+5uFIDRo0fz+OOPs3XrVmbPnp1j8am8fn4FfTaye+ihh5g/fz7z5s3jgQceuOTPSFUv+mzt2bOHl19+mc8++4wffviBPn362OJYFVhcXBytWrVi27YiTYwyVYRdAciloG/sNWvWLLC+fv36hfrGn92pU6cICAigZs2aJCYmsmHDBgCuvfZaxowZw7Fjx6hTpw7vv/8+7du3z7HtFVdcQUBAAOvWreOGG27gnXfecV8NKIzatWtz+vTpPOt2795NeHg44eHhrF+/nsTERG699VYmTpzIsGHD8Pf358CBAwWuBJeZmcn+/fu58cYbuf7661m4cCFnzpy5KIb8FtXo0aMHDz74IHPnzmXIkCF5xp972xEjRtCpUycaNWpEaGgoAAcOHODee+/ls88+y/HegwcP8rvf/Y577rkHf39/5s+fz/jx4zl69Cjr16+nS5cupKWlsXPnTve+LsfkyZOZPPm3ySmnTp0iMDAQgLfeeuuS2+f32cjt9ttvZ9KkSaSlpblH/hf0c+3ZsyfTp09n2rRpAJw4cYLTp09Tq1YtrrjiCg4fPszHH3+c59UwU/4lJCQwePBgWrRoYcufmwJZAuBhvXr1YtasWURERBASEkLnzp0BaNy4MS+88AJdunShcePGdOzYkYyMDPd2Wd/g3nrrLR555BHOnTtHy5YtmTdvXqGP3bdvXwYPHsyyZct4/fXXueGGG9x106ZN44svvsDb25t27drRu3dvfH192bFjB126dAGcg93effddvL2989x/RkYG99xzD6dOnUJVefLJJ92rcWUZOnQoDz/8MLGxse5Bi1m8vb2JiYlh/vz5eTaY9erVo2vXroSFhdG7d2+mTp1Kw4YNadu2bY7uhqSkpDy7ArZu3crTTz+Nl5cX1apVY+bMmVSvXp0lS5bwxBNPcOrUKdLT0/njH/9YpAQgtxdeeIE77riDwMBAOnfuzJ49ewp8f36fjdyqV6/OjTfeSN26dd3/FxEREfj4+NC+fXtGjBjh7rYBeO655xg1ahRhYWF4e3vz/PPPM3DgQCIjIwkNDaVly5Z07dq12Odryl5SUhJ9+vShdu3arFy5slBXzUzVVWFXAiyK/FYCrGirmY0ePZqOHTte8hJ3VXTu3DnCw8PZtGkTV1xxBQDTp0+nWbNm7v77yiYzM5OOHTvy/vvvu2/GUhYq4u9OZXb27Fl69OjB9u3bWbduHR07dvR0SKYcqJQrAVZVEydO5Jtvvqm0jVlxrF69mjZt2jB69Gh34w/w+OOPV9qf1/bt27n66qu56aabyrTxN+VPZmYmjRs3ZtGiRdb4m0KxKwD2LcaYIrHfnfIjNTWV6tWr5znA01RtdgXAGGMqqZkzZ3Lddddx/Phxa/zNZbEEAOe0KGNM4dnvTPnwySefMHr0aBo1apSj28uYwqjyCYCfnx/Hjh2zP2jGFJKqcuzYMfz8/DwdSpX2ww8/cOeddxIeHs6iRYvynY1jTH48Pg1QRO4AXgDaAp1UdWO2ugnAg0AG8ISqfuoqjwLmAzWAVcAYLWILHhQUhMPh4OjRo8U5DWOqFD8/P4KCgjwdRpWVlJRETEwMderUYcWKFfnef8KYgng8AQASgIHA7OyFItIOGAqEAk2A1SLSWlUzgJnASGADzgSgF/BxUQ5erVq1HEuyGmNMeXf+/HkaNmzI3Llz3YtLGXO5PJ4AqOoOyHPZ2/7AIlW9AOwRkV1AJxHZC9RR1fWu7d4GBlDEBMAYYyqKjIwMvLy8aNGiBd9++60N+jPFUp7HAAQC+7O9drjKAl3Pc5cbY0yl9vTTTzNs2DAyMjKs8TfFViYJgIisFpGEPB79C9osjzItoDy/Y48UkY0istH6+Y0xFdWMGTN49dVXadCggQ34MyWiTLoAVPXmImzmAJpmex0EHHSVB+VRnt+x5wBzwLkQUBHiMMYYj1q1ahWjR4+mb9++vPLKK54Ox1QS5bkLYDkwVER8RaQF0Ar4VlWTgBQR6SzOa2D3Ass8GagxxpSW77//niFDhtC+fXsWLlxo3/5NifF4AiAit4uIA+gCrBSRTwFUdRuwGNgOfAKMcs0AAHgUeBPYBezGBgAaYyqpU6dOERwczEcffWTT/UyJqvL3AjDGmPIoMzMTLy+vi54bcznsXgDGGFOBZGRkMGjQIKZOnQpgjb8pFfapMsaYcuapp55i6dKl1KxZ09OhmErMEgBjjClHpk+fzrRp0xgzZgyjRo3ydDimErMEwBhjyomVK1cyZswY+vbty7/+9S9Ph2MqOUsAjDGmnHA4HERFRdl0P1MmbBaAMcZ4mKq6l/ZNT0/Hx8fjt2kxlYTNAjDGmHLqzJkzdO/enVWrVgFY42/KjCUAxhjjIRkZGdx11118/fXXdnMfU+Ys1TTGGA8ZO3YsK1asYMaMGfTu3dvT4Zgq5pJXAETERqIYY0wJi42NJTY2lieffJJHH33U0+GYKqgwXQBvlnoUxhhTxezYsYP+/fu7V/szpqwVpgvAW0T+qqqTsheKSG3g/6nqHaUTmjHGVD5ZI/5nzJhBWlqaTfczHlOYKwAPAF1E5KGsAhFpD2wEdpRWYMYYU9k4HA66d+/Ojh07EBGqV6/u6ZBMFXbJBEBV04GBwKMicpuIPAwsw3l73kkFb31pIjJVRBJF5AcR+Y+I1M1WN0FEdonIjyJya7byKBHZ6qqLFRs+a4wp51JSUujbty9btmwhLS3N0+EYU6hBgDOAYcA/gXnAEOBaVV1dQjH8FwhT1QhgJzDBddx2wFAgFOgFzMg2IHEmMBJo5Xr0KqFYjDGmxKWnp3PXXXexdetWFi9eTEREhKdDMqZQXQBbgDDgMaAa0BqYKyKTRWRocQNQ1TjXVQaADUCQ63l/YJGqXlDVPcAuoJOINAbqqOp6dS5j+DYwoLhxGGNMaRk7diwrV67k9ddfp1cv+75iyodLDgJU1TnZX4tIEBABhAN9gEUlGM8DwP+5ngfiTAiyOFxlaa7nucuNMabcOX/+PJs3b+ZPf/qTTfcz5cplLwSkqg6cje6qwm4jIquBRnlUPauqy1zveRZIBxZkbZbX4Qsoz+/YI3F2F9CsWbPChmyMMSXCz8+P1atX2xK/ptwp1FLAIuIrIhNFZJWIvO8anHdVYQ+iqjeralgej6zG/z4gBhimv92dyAE0zbabIOCgqzwoj/L8jj1HVaNVNbpBgwaFDdkYY4pl06ZNxMTEcOLECXx9fW26nyl3CnsvgKlAc5wD9BxAPWCeiNxT3ABEpBcwDuinqueyVS0HhrqSjxY4B/t9q6pJQIqIdHaN/r8X56wEY4wpFxwOB3379mXr1q1cuHDB0+EYk6fCXpNqD/RX1ZMi0klVu4qIL7AGeLeYMUwHfIH/umbzbVDVR1R1m4gsBrbj7BoYpaoZrm0eBeYDNYCPXQ9jjPG4lJQUYmJiOHPmDF9//TWNGuXV+2mM5xU2AZjKb1cLkkXkFSARyCxuAKp6dQF1k4HJeZRvxDkzwRhjyo309HSGDh1KQkICq1atIizM/kyZ8qsw6wB0AVaq6nFX0R04VwH0dz03xhgDJCUlsXXrVt544w169uzp6XCMKVBhrgDcB7whIjuBT4BPVHVh6YZljDEVT9OmTdm2bRu1a9f2dCjGXFJhlgJ+RFU7Ai8AAcB8EVkvIv8QkW52u2BjTFW3fPlyHn/8cdLT063xNxVGYWcBoKqJqvqqqvYCfg98hbML4JvSCs4YY8q7+Ph47rrrLr799ltSU1M9HY4xhVbolSlE5Pc47wlwEkgAfgCeUlWb42KMqZL2799P3759qV+/PsuXL6dmzZqeDsmYQrucpaneBUa5tonAuf5+KJDvKH5jjKmsTp8+TZ8+fTh79qxN9zMV0uUkALtU9T+u5++XRjDGGFNR/PDDD+zbt48lS5bYdD9TIV0yARCRt4FNwHoR+ZOq/qv0wzLGmPLt+uuvZ+/evQQEBHg6FGOKpDCDAN9yva8RMFxE9onIchH5m4jYOgDGmCrltddeY84c501SrfE3FVlhpgF+pqqvqOp9qtoBuAp4DtgJXFvK8RljTLmxdOlSnnzySeLi4vjtvmXGVEyFGgPgWvf/GaALcBZnl8BiVX2nFGMzxhiPO336NFOmTGHt2rVs2LCBTp068c477+C6d4kxFVZx7wY4rLQCM8aYsnbs2DGWLVvGn/70J2JjYwHw8/Nj+vTppKWlMXbsWFasWEGNGjU8HKkxxVcSdwNcUGrRGWNMGZg4cSJLly4lISEBAF9fX+6//34AqlevzpEjR/D19fVkiMaUOI/fDdAYY8qCqrJv3z7Wrl3L2rVr+eWXX4iLiwOcC/o0adKEoUOH0q1bNzp16pSjwbfG31RGhUoAVHVFtpd3AINxzgoo9iwAEfkb0B9nMnEEGKGqB111E4AHgQzgCVX91FUeBcwHagCrgDFqI3KMMdlk/UkQEebOncvf/vY39u/fDzhH73fr1o0LFy7g6+vL/PnzPRipMZ5R6HsBZFHVVFVd6JoZcLAEYpiqqhGuGQYrgEkAItIOGIpztcFewIxsNx6aCYwEWrkevUogDmNMBZaRkcGWLVuIjY1l0KBBNGzYkMTERADq1atHly5dmD59Oj/88APJycksXbrUvtmbKu1yVgIsFap6OtvLWkDWN/n+wCLXvQb2iMguoJOI7AXqqOp6cC9UNAD4uMyCNsZ4XFpaGhcuXMDf358NGzbQq1cvTp06BUBwcDC9e/d2v3fgwIEMHDjQU6EaUy55PAEAEJHJwL3AKeBGV3EgsCHb2xyusjTX89zl+e17JM6rBTRr1qzkgjbGlKlff/2Vb775xt2Hv379eiZNmsS4ceNo1aoVQ4YMoVu3btxwww32u25MIZRJAiAiq3GOGcjtWVVdpqrPAs+6+vwfB54H8ppkqwWU50lV5wBzAKKjo22cgDEVxOnTpzl8+DCtWrUiLS2Nhg0bkpKSgogQERHBQw89RNeuXQHnJf7Zs2d7OGJjKpYySQBU9eZCvnUhsBJnAuAAmmarCwIOusqD8ig3xlRgx44d46uvvuLLL79k7dq1bN68mc6dO/P1119TrVo1XnzxRZo3b07Xrl1tCV5jSsBlDwIsaSLSKtvLfjinFwIsB4aKiK+ItMA52O9bVU0CUkSksziX4roXWFamQRtjCmfBAggOBi8v578Lfls25ODBg3z00Ufu1yNGjGDAgAHMmDEDf39/nnvuOf7+97+760eNGkVMTIw1/saUEPH07DkR+QAIwTkNcB/wiKoecNU9CzwApAN/VNWPXeXR/DYN8GNgdGGmAUZHR+vGjRtL4zSMMbktWAAjR8K5c+6i9OrV+Xfnzkw9cIDdu3cDcOTIERo0aMCGDRtIT0/nmmuusdH5xpQQEYlX1eg86zydAJQlSwCMKRsnT57ENySEGkeOXFS338uL0X370q1bN7p160aHDh3w8SkX45GNqXQKSgDst84YUyyqiojw888/M378eOLj4/n555/JyOf9QaosXbq0LEM0xuTB42MAjDEVR3JyMp9++ikvvvgigwcPpkWLFvzrX/8CoEaNGnz33XdERkbyj3/8gwtXXpnnPsSm6BlTLtgVAGNMng4fPsymTZvw8fHhlltuITU1lSZNmpCWlgbAVVddRadOnQgJCQGgcePG7Nmz57cdNGt20RgAataEyZPL8jSMMfmwBMAY4xYbG8tnn31GfHw8Bw4cAKBbt27ccsstVK9enTlz5tCsWTM6duxI3bp1C97ZMNfdwp99Fn75xZkQTJ78W7kxxqNsEKAxVczBgweJj493P1JSUlizZg0AMTEx7N69m6ioKDp27EhUVBSRkZHUqVPHs0EbY4rEBgEaUwWpKg6Hg82bN9O3b19EhDFjxhAbGwuAl5cXbdq0ITo6mszMTLy8vFi2bBne3t6X2LMxpjKwBMCYSiQhIYH33nuP+Ph4Nm3axNGjRwHYtWsXV111Ff369ePqq6+mY8eOdOjQgVq1auXY3hp/Y6oOSwCMqWBUlb1797ov4W/atIm///3vXHPNNfz4449MmTKF0NBQYmJi3Jfxg4Kcq2ffdNNN3HTTTR4+A2NMeWAJgDHlmKry888/4+vrS1BQEAkJCXTr1o0TJ04A4OPjQ1hYmPs2uH369CElJYUaNWp4MmxjTAVgCYAx5Uh6ejpz584lMTGRrVu3smnTJk6dOsWECRP4xz/+QXBwMIMHD3YP0gsPD8fPz8+9ffbnxhhTEJsFYEwZ27ZtG9u3bycxMZEff/yRxMREOnXqxIwZM1BV6tatS0ZGBm3btiUqKoqoqCi6d+9O69atPR26MaaCsVkAxpQhVeXQoUM5GnhfX1+mTJkCwJ133sn27dsBaNasGSEhIVx99dUAiAg//fQTDRo0wHmzS2OMKR2WABhTROfPn2fXrl0kJiZy+PBhRo0aBcDtt9/OsmW/3aG6Zs2adOvWzf169uzZ+Pv706pVq4tG4QNcmc8SusYYU5LKTReAiDwFTAUaqGqyq2wC8CCQATyhqp+6yqP47XbAq4AxdjtgUxpUlcOHD/Pjjz9yww034OXlxbRp04iNjWXv3r1kfex8fHw4e/Ys1atXZ8mSJRw6dIiQkBDatGlDYGAgXl522w1jTNkr910AItIUuAX4JVtZO2AoEAo0AVaLSGtVzQBmAiOBDTgTgF7Ax2Udt6k8Lly4gLe3Nz4+Pqxbt44333zTfQk/a4T9vn37aNasGfXq1aNTp07ce++97ka+devWVK9eHYDBgwd78lSMMaZQykUCALwKPAMsy1bWH1ikqheAPSKyC+gkInuBOqq6HkBE3gYGYAmAKaSkpCRWrFhBYmKiu5Hfs2cP69at47rrrsPhcPD5558TEhLCPffc427k69evD8Dw4cMZPny4h8/CGGOKx+MJgIj0Aw6o6ve5Bj0F4vyGn8XhKktzPc9dbgzgnEq3c+dO9wC8rH8nTJhA//792b17NyNHjsTPz4/WrVsTFRXF3XffTcOGDQEYOnQod911l4fPwhhjSleZJAAishpolEfVs8CfgZ55bZZHmRZQnt+xR+LsLqCZ3Ye8UsjIyGDr1q04HA4cDgf79+/H4XDQu3dvhg4disPhIDQ01P3+Jk2a0KZNG6pVqwZAdHQ0e/fupWnTpnn2zdvoe2NMVVAmCYCq3pxXuYiEAy2ArG//QcAmEemE85t902xvDwIOusqD8ijP79hzgDngHARY9LMwZWnlypXs27fP3bg7HA5uvPFGJk2aRGZmJlFRUWRmZgLO9eubNGlCeHg44Ez03n33XUJCQggJCaF27do59u3n50fz5s3L/JyMMaY88WgXgKpuBdxznlz9+9Gqmiwiy4GFIvIKzkGArYBvVTVDRFJEpDPwDXAv8HrZR28uR2ZmJikpKVxxxRUAvP3222zbts3duDscDqKioli8eDEAf/jDHzhw4AA+Pj4EBgYSFBREzZo1AahWrRpLly7lyiuvJCgoiEaNGuW4iY2XlxfD7J7zxhhTII+PAciPqm4TkcXAdiAdGOWaAQDwKL9NA/wYGwDoURkZGRw+fJjjx48TFhYGwBtvvMFXX33lbtwPHDhAaGgomzdvBmDGjBls3ryZoKAggoKC6Ny5M507d3bvMy4ujoCAABo2bJjnZfq+ffuWzckZY0wlVa4SAFUNzvV6MjA5j/dtBMLKKKwqLT09nUOHDuFwODh06BADBgwA4NVXX2XJkiU4HA4OHjxIeno69erVIzk5GYANGzawceNGgoKCuP766wkKCsqxlG1cXBz+/v75zo9v165dqZ+bMcZUZeUqATBlJzU1leTkZE6cOMHx48dJSkpi//79PPbYY9SoUYNp06bx8ssvk5SU5O5rBzh37hw1atTgwoUL+Pn50aNHD/e3+KCgIFQVEeGdd94p8Ph16tQp7VM0xhhTgHKzEmBZqEwrAaoqZ8+e5cSJEzke119/PfXr12fjxo289dZbF9UvXbqUkJAQYmNjGTNmzEX73blzJ61ateKDDz5gxYoVNG3aNEcDHxoamqO/3RhjTPlV7lcCrKpUlZSUlBwNdOvWrQkMDOSXX35h9uzZnDx5Mkf9Sy+9RI8ePVi5cmWe/eD//e9/ufnmm9m3bx8LFiwgICCAunXrEhAQQGBgoLvx/v3vf8+sWbMICAggICCAJk2aEBQU5P5mPmjQIAYNGlSmPw9jjDFlxxKAYsrMzOT06dPuBrpevXo0b96clJQUZs6cedE38Iceeog777yTxMREwsLCyMjIyLG/2bNnM3LkSI4ePcqUKVPcjXdWQ57VZx4WFsaUKVPcdVmPrH72SzXgYWFh7gF7xhhjqh5LAIro/PnzNGnShJMnT5K9G2X8+PG8+OKLpKWlMW7cOHx8fHI00KmpqQA0bNiQcePGXdSAZw1+i4yMJC0tLd9FaYKDg3nmmWdK/0SNMcZUSpYAFJGvry/Dhw/niiuuyPEtPasBDwgIICUlhVq1auXZiAcEBDB58kUTHNzs7nHGGGNKkyUARSQivPbaawXW+/v7l2FExhhjTOHZ10xjjDGmCrIEwBhjjKmCLAEwxhhjqqAqtRCQiBwF9pXgLusDySW4P0+qLOdSWc4D7FzKq8pyLpXlPMDOpSDNVbVBXhVVKgEoaSKyMb8VliqaynIuleU8wM6lvKos51JZzgPsXIrKugCMMcaYKsgSAGOMMaYKsgSgeOZ4OoASVFnOpbKcB9i5lFeV5Vwqy3mAnUuR2BgAY4wxpgqyKwDGGGNMFWQJQBGIyP8TkSMikuDpWIpDRJqKyBciskNEtonIGE/HVFQi4ici34rI965z+YunYyoOEfEWkc0issLTsRSXiOwVka0iskVENno6nqISkboiskREEl2/M108HVNRiEiI6/8i63FaRP7o6biKSkSedP3OJ4jIeyLi5+mYikJExrjOYVtZ/X9YF0ARiEg34AzwtqpW2HvqikhjoLGqbhKR2kA8MEBVt3s4tMsmzjsu1VLVMyJSDfgKGKOqGzwcWpGIyFggGqijqjGejqc4RGQvEK2qFXqetoi8BaxT1TdFpDpQU1VPejisYhERb+AAcK2qluQaKWVCRAJx/q63U9VfRWQxsEpV53s2sssjImHAIqATkAp8Ajyqqj+V5nHtCkARqOpa4Lin4yguVU1S1U2u5ynADiDQs1EVjTqdcb2s5npUyOxWRIKAPsCbno7FOIlIHaAb8G8AVU2t6I2/y03A7orY+GfjA9QQER+gJnDQw/EURVtgg6qeU9V04Evg9tI+qCUABgARCQYigW88HEqRuS6bbwGOAP9V1Yp6LtOAZ4BMD8dRUhSIE5F4ERnp6WCKqCVwFJjn6pp5U0RqeTqoEjAUeM/TQRSVqh4AXgZ+AZKAU6oa59moiiQB6CYi9USkJnAb0LS0D2oJgEFE/IEPgD+q6mlPx1NUqpqhqh2AIKCT67JahSIiMcARVY33dCwlqKuqdgR6A6NcXWgVjQ/QEZipqpHAWWC8Z0MqHlc3Rj/gfU/HUlQiEgD0B1oATYBaInKPZ6O6fKq6A5gC/Bfn5f/vgfTSPq4lAFWcq7/8A2CBqn7o6XhKguvS7Bqgl2cjKZKuQD9Xv/ki4Pci8q5nQyoeVT3o+vcI8B+c/ZwVjQNwZLuqtARnQlCR9QY2qephTwdSDDcDe1T1qKqmAR8C13k4piJR1X+rakdV7Yazi7lU+//BEoAqzTVw7t/ADlV9xdPxFIeINBCRuq7nNXD+YUj0aFBFoKoTVDVIVYNxXp79XFUr3DeaLCJSyzXAFNcl8544L3dWKKp6CNgvIiGuopuACjdYNpe7qMCX/11+ATqLSE3X37ObcI5lqnBE5ErXv82AgZTB/41PaR+gMhKR94AeQH0RcQDPq+q/PRtVkXQFhgNbXX3nAH9W1VWeC6nIGgNvuUY1ewGLVbXCT6GrBBoC/3H+bcYHWKiqn3g2pCIbDSxwXTr/Gbjfw/EUmauf+RbgD56OpThU9RsRWQJswnnJfDMVd1XAD0SkHpAGjFLVE6V9QJsGaIwxxlRB1gVgjDHGVEGWABhjjDFVkCUAxhhjTBVkCYAxxhhTBVkCYIwxxlRBlgAYY4wxVZAlAMYYY0wVZAmAMabEiEgNEfnSdWOmYBEp0qp/IlJdRNa67vBmjCkFlgAYY0rSA8CHqppRnJ2oairwGTCkRKIyxlzEEgBjTKGIyBcicovr+d9FJDaPtw0DluWxbUvXbXSvcV0ZSHTdUjdBRBaIyM0i8rWI/CQiWTcLWuranzGmFFgCYIwprOeBZ0VkGBAJPJm90rVGfktV3ZurPATnHSfvV9XvXMVXA68BEUAb4G7geuAp4M+u9yQA15TKmRhj7GZAxpjCUdW1rjuujQV65HGZvz5wMldZA5xXBAap6rZs5XtUdSuAiGwDPlNVFZGtQLDreBkikioitVU1peTPyJiqza4AGGMKRUTCcd518UI+DfKvgF+uslPAfpx3nszuQrbnmdleZ5Lzi4kvcL6oMRtj8mcJgDHmkkSkMbAA6A+cFZFbc7/HdftSbxHJngSkAgOAe0Xk7ss8Zj3gqKqmFTlwY0y+LAEwxhTIde/4D4E/qeoO4G/AC/m8PQ5nX76bqp4FYoAnRaT/ZRz6RmDVZQdsjCkUUVVPx2CMqSREJBIYq6rDS2BfHwITVPXH4kdmjMnNrgAYY0qMqm4GvhAR7+LsxzWjYKk1/saUHrsCYIwxxlRBdgXAGGOMqYIsATDGGGOqIEsAjDHGmCrIEgBjjDGmCrIEwBhjjKmCLAEwxhhjqiBLAIwxxpgq6P8DvUtb9fyaKdoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dhdK = np.empty(0)\n",
    "for xp in X[::1000]:\n",
    "    M0 = timml.Model3D(kaq=K1, z=[0., -b])\n",
    "    timml.inhomogeneity1d.StripInhom3D(M0, x1=2*L/3., x2=L, kaq=K2, z=[0,-b], topboundary='conf')\n",
    "    timml.ImpLineDoublet1D(M0, xld=0.)\n",
    "    timml.HeadLineSink1D(M0, xls=L, hls=0.)\n",
    "    timml.LineSink1D(M0, xls=xp, sigls=-1.)\n",
    "    M0.solve(silent=True)\n",
    "    A = M0.headalongline(X, 0.).flatten() \n",
    "    dhdK = np.append(dhdK, np.sum(A*b*-R/(K1*b)))\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "plt.subplots(figsize=[8,2])\n",
    "plt.plot(X[::1000][1:]/1000., dhdK[1:], 'k--', mfc='none', label='adjoint sensitivity, semi-analytical')\n",
    "plt.plot(ocol/1000., dhdK[int(ocol/1000)], 'ro')\n",
    "plt.xlabel('$x$ (km)')\n",
    "plt.ylabel(r'$\\partial h / \\partial K$')\n",
    "plt.legend(loc=2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src onerror=\"var cell = this.closest('.code_cell');var editor = cell.querySelector('.input_area');editor.style.background='rgba(0, 200, 0, 0.2)';this.parentNode.removeChild(this)\">"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-325.020001 \n"
     ]
    }
   ],
   "source": [
    "set_background('rgba(0, 200, 0, 0.2)')\n",
    "print('%.6f '% (dhdK[int(np.where(X[::1000]==float(ocol))[0])]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3c. Numerical, continuous"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfsAAACdCAYAAACzUdEhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlAUlEQVR4nO3deVyVZfr48c8limSIG2UoCFaWK+5Lai6NpuMoaqk1P8alza/TmJXTTFrfUbOxfPk1rbExv5Vru5UaNv7KMmlxydxyxcwEJTEUFdxAwev3xzmcHwdQQZaHA9f79eLFOc/9LNd5WK7nvp/7uW9RVYwxxhhTflVyOgBjjDHGlCxL9sYYY0w5Z8neGGOMKecs2RtjjDHlnCV7Y4wxppyzZG+MMcaUc5WdDqCkBAcHa0REhNNhGGOMMaVmy5Ytx1X1htzLHUn2IhIGLAFuAi4Br6vqKyJSG/gAiADigWGqetK9zUTgISALGKeqn1/pGBEREWzevLnEPoMxxhhT1ohIQn7LnWrGzwT+qqpNgE7AX0SkKTABWKOqjYA17ve4y+4HmgF9gbki4udI5MYYY4yPcSTZq2qSqm51vz4N7AXqAwOBxe7VFgOD3K8HAu+raoaqHgR+BjqUatDGGGOMj3L8nr2IRACtge+BuqqaBK4LAhG50b1afWBjjs0S3ctKzZw5c1i5cqXXMn9/fz799FMAZsyYwZo1a7zKa9SowdKlSwGYOnUq69ev9yq/6aabWLRoEQATJ05k+/btXuUNGzZk7ty5ADz55JPs27fPq7xp06bMnDkTgD//+c8cOnTIq7xt27ZMnToVgAceeIBjx455lXft2pUJEyYAcP/993P27Fmv8l69evH4448DMGjQIC5duuRVPmDAAB555BEyMjIYNmwYuQ0dOpQ//elPpKamMmLEiDzlI0aM4N577+Xo0aOMGTMmT/no0aPp168fBw8eZPz48XnKx40bR8+ePdm7dy/PPvtsnvKnn36ajh07snXrVqZNm5anfPLkyURGRrJ+/XpmzZqVp/zFF1+kUaNGfPXVV7z22mt5ymfNmkVYWBirVq3y/Bxzmjt3LsHBwSxbtowPPvggT/n8+fMJDAzk3XffJSYmJk/5O++8g5+fHwsWLOCLL77wKvP392fx4sWe43z33Xde5UFBQcybNw+A2bNn57mldeONNzJ79mzP59y9e7dXeYMGDXjhhRcAmDJlCgcOHPAqv+222/jHP/4BwIQJEzhy5IinLDg4mClTphAUFJTnMxljnOFosheRQOBj4AlVTRORy66az7I8g/qLyGhgNLj+WRWn9PR0zpw547XM39/fqzwtLc2r3M/v/99pOHv2LKdOnfIqr1atmuf16dOnOXHihFd57dq1Pa9PnTrF8ePHvcpTU1M9r1NSUkhOTvYqz3m85ORkfvvtt8uWHz16lNOnT1+2/NdffyUrKyvfclXNc6GRM76srKx8y7PPV2ZmJvHx8XnKs+O5cOECv/zyS57y7J/H+fPn+fnnny9bfvbsWX766ac85dkXN2lpacTFxeUpP3/+POD6nHv27MlTnpGRAbjO/a5du/KUX7hwAYDffvuNHTt25CnPPp9HjhzJc6EHrvMKcOjQIbZu3epVVrVqVc/r+Pj4PMk85+/O/v372bRpk1d5aGio5/XevXvZuHGjV/nJkyc9r3fs2MGPP/7oVZ792QG2bt3qdf4PHTrEhQsXePXVV/N8JmOMM8SpiXBEpArwKfC5qs5yL9sH9HDX6kOAWFW93d05D1V90b3e58AUVd1wuf23a9dOrYOeMaVv5cqV3HnnndSsWdPpUIypcERki6q2y73ckXv24qrCzwf2Zid6txhgpPv1SOCTHMvvF5GqItIQaAR4V1WMMWXCgAEDqFmzJunp6XluLxhjnOFUb/wuwHDgLhHZ7v7qB0wHeovIfqC3+z2quhtYCuwBPgP+oqpZ+e/aGFMWPP300/Tq1Yvvv//e6VCMqfAca8YvadaMb4yzjh8/TocOHUhPT+eHH36gfv1S7VNrTIVUpprxjTHlX3BwMDExMZw+fZpBgwZ5OjwaY0qfJXtjTIlp3rw577zzDlu2bOGJJ55wOhxjKizHn7M3xpRvUVFR/O///i/dunVzOhRjKixL9saYEvfII48ArrEDDhw4wK233upwRMZULNaMb4wpNS+88AKtW7fOM2KfMaZkWbI3xpSakSNHEhgYSFRUFCkpKU6HY0yFYcneGFNqQkNDWb58Ob/++itDhw7l4sWLTodkTIVgyd4YU6o6derEG2+8wdq1a/Od4MgYU/ysg54xptQNHz6cgwcP0qGDzVRtTGmwZG+MccSkSZM8r9PS0mxKXGNKkDXjG2MctXjxYho1apTvNMbGmOJhyd4Y46jOnTtz8eJFBg4cyOnTp50Ox5hyyZK9McZRjRo1YunSpezdu5fhw4dz6dIlp0MyptyxZG+McVyvXr2YPXs2n3zyide9fGNM8bAOesaYMmHs2LHs3r2bWrVqOR2KMeWOJXtjTJkgIrz22muICABZWVn4+fk5HJUx5YMjzfgiskBEkkVkV45lU0TkVxHZ7v7ql6Nsooj8LCL7RKSPEzEbY0pedqJfs2YNLVu2JCkpyeGIjCkfnLpnvwjom8/y2arayv21CkBEmgL3A83c28wVEbvcN6YcCw4OJj4+nsGDB5Oenu50OMb4PEeSvap+A5wo4OoDgfdVNUNVDwI/AzbsljHlWMuWLXnrrbf4/vvvGT16NKrqdEjG+LSy1ht/rIjscDfzZ/fSqQ8czrFOonuZMaYcGzx4MFOnTuWtt97ipZdecjocY3xaWUr2rwG3AK2AJCD7r1vyWTffy3wRGS0im0Vk87Fjx0okSGNM6fnv//5vhg4dyp49e6x2b0wRlJne+Kr6W/ZrEXkD+NT9NhEIy7FqKHDkMvt4HXgdoF27dvafwRgfJyK8/fbbVKlSxdN5zxhTeGWmZi8iITneDgaye+rHAPeLSFURaQg0AjaVdnzGGGf4+/sjIuzbt4/+/ftz4kRBu/sYY7JdU7IXkeuL0iNeRN4DNgC3i0iiiDwEzBCRnSKyA+gJPAmgqruBpcAe4DPgL6qada3HNsb4ppSUFFavXs19991HZmam0+EY41OkIPfBRKQSrsffooH2QAZQFTgGrAJeV9X9JRhnobVr1043b97sdBjGmGK0cOFCHnzwQcaNG8crr7zidDjGlDkiskVV2+VeXtB79muBL4GJwC5VveTeaW1ctfDpIrJcVd8uroCNMSa3Bx54gJ07dzJ79mxatGjBww8/7HRIxviEgtbsq6jqxaKuU5qsZm9M+ZSZmUn//v05d+4csbGxVKpUZroeGeO4ItXsC5LEy1KiN8aUX5UrV2bp0qX4+/tbojemgK76lyIivUXkDRFp5X4/usSjMsaYKwgKCiIgIIBTp07x+OOPc+bMGadDMqZMK8hl8aPA34A/ichduAa9McYYx23bto1XX32VESNGcOnSJafDMabMKkiyP6aqp1T1KeBuXL3xjTHGcT179mTmzJksX76c5557zulwjCmzCpLs/5P9QlUnAEtKLhxjjCmcJ554glGjRjF16lQ+/PBDp8MxpkwqUG98X2S98Y2pODIyMrjrrrs4duwYu3fvpkqVKk6HZIwjivqc/ZV2vE5VuxR1P8YYc62qVq3KsmXLyMrKskRvTD6K47mVesWwD2OMKZK6detSr149srKymDNnDhkZGU6HZEyZUaBkLyJz3NPH3iEi1XMVl8/7AMYYn/Ttt98ybtw4xowZY9PiGuNW0Jr9TiASmA7Ei8hBEYkRkWlA7uRvjDGO6dGjB5MmTWLRokW8/PLLTodjTJlQ0BH0Xs/5XkRCcSX/FsDnJRCXMcZcs8mTJ7Nr1y6eeuopmjZtSp8+fZwOyRhHFXRs/DuAjepDbWLWG9+Yiu3MmTN06dKF5ORkfvnlF6677jqnQzKmxBW1N/5IYK6I7MM1p/xnqnq0OAM0xpjiFBgYSExMDMnJyZboTYVXoHv2qjpGVVsDU4BawCIR2SAiL4hINxHxK8kgjTHmWoSHh9O+vWvQz88//5zMzEyHIzLGGQXtjV8ZQFXjVHW2qvYF7gK+A4YC3xfmoCKyQESSRWRXjmW1ReQLEdnv/l4rR9lEEflZRPaJiN18M8YUyqZNm+jbty9///vfnQ7FGEcUtDf+ptwLVPW8qq5S1cfyuz9wFYuAvrmWTQDWqGojYI37PSLSFLgfaObeZq61JBhjCqNDhw489thjzJ49m4ULFzodjjGlrqDJXorzoKr6DXAi1+KBwGL368XAoBzL31fVDFU9CPwMdCjOeIwx5d+sWbPo1asXY8aMYf369U6HY0ypKmgHvRtEZPzlClV1VjHEUldVk9z7SxKRG93L6wMbc6yX6F6Wh4iMBkYDNGjQoBhCMsaUF5UrV+aDDz6gY8eO3HPPPfz0008EBQU5HZYxpaKgyd4PCKSYa/gFlN8x830E0D0ewOvgevSuJIMyxvie2rVrs3LlSrZu3WqJ3lQoBU32Sao6tUQjgd9EJMRdqw8Bkt3LE4GwHOuFAkdKOBZjTDnVuHFjGjduDEBcXBy33XYblSoVxzQhxpRdjtyzv4wYXM/z4/7+SY7l94tIVRFpCDQinw6DxhhTGLt376ZVq1b885//dDoUY0pcQZP9XVdbQUQKfEEgIu8BG4DbRSRRRB7CNe5+bxHZD/R2v0dVdwNLgT24BvT5i6pmFfRYxhiTn6ZNmzJs2DAmT57Mxx9/7HQ4xpSogg6XGwt8DHyiqodyLPcHuuKqia9V1UUlE2bh2XC5xpirSU9Pp2fPnuzYsYP169fTsmVLp0MypkguN1xuQWv2fYEs4D0ROSIie0TkF2A/8EdgdllK9MYYUxABAQEsW7aMWrVqERUVRUpKitMhGVMiCjrrXTowF9eANlWAYOC8qp4qwdiMMabEhYSEsGLFCpYtW0bNmjWdDseYElHQ3vgeqnoRSCqBWIwxxhHt2rWjXTtXy+fJkyepWbMmheiGZEyZZ8+bGGOMW0JCAk2bNmXOnDlOh2JMsbJkb4wxbmFhYXTs2JEnn3ySL774wulwjCk2luyNMcatUqVKvPXWW57H8vbv3+90SMYUi6smexG53v09sOTDMcYYZ1WvXp2YmBj8/PyIiooiNTXV6ZCMKbKC1OxrichYXM/TG2NMudewYUM+/vhjbr/9dqdDMaZYFCTZ/w4YBdycYya6fInI08URlDHGOK179+6sWLGCGjVqkJmZ6XQ4xhRJQZL9JuBBIEFVk3MWiMjSHF8fAg+XRJDGGOOU48eP07lzZ5YsWeJ0KMZcs6s+Z6+qe90vd+RTnKaqngQvIq8VV2DGGFMW1KhRg8DAQB555BFuu+02OnXq5HRIxhRaoXrji8h7IvK0iPxeRG4Ack8X9WzxhWaMMc6rUqUKH374IaGhoQwaNIjExESnQzKm0Ar76N084BxwL/AFMFtEqmcXquqJYozNGGPKhDp16hATE8PZs2cZNGgQ586dczokYwqlsMm+AfA1MEZVWwEfAc8Vd1DGGFPWNGvWjHfffZfz589z/Phxp8MxplAKOzb+bcBQoKmInMZ1H7+7iPwH2KGqx4o7QGOMKSsGDBhA3759qVKlitOhGFMoha3Z/0tVo1T1VuAu4G2gGhANrCqOgEQkXkR2ish2EdnsXlZbRL4Qkf3u77WK41jGGFNYVapU4fz584wYMYJPPvnE6XCMKZDCJvvV7mT8f4HJwF+B91T1QVVtX4xx9VTVVqrazv1+ArBGVRsBa9zvjTHGMXFxcURHR7Nz506nQzHmqgqV7FW1NXAL8BTwDfAqroRf0gYCi92vFwODSuGYxhiTr+uuu47ly5cTFBREVFSU3cM3ZV6hJ8JR1SxV3a2qy1T1U1Ut7qGlFFcLwhYRGe1eVldVk9zHTwKuOJKfMcaUtPr167NixQqSkpIYMmQIFy5ccDokYy6rLM5610VV2wC/B/4iIt0KuqGIjBaRzSKy+dgx6ytojClZHTp0YP78+fz444/s27fP6XCMuSxRVadjuCwRmQKcAR4BeqhqkoiEALGqesUZKtq1a6ebN28uhSiNMRVdSkoKderUcToMYxCRLTn6u3mUqZq9iFyfPUiPe2rdu4FdQAww0r3aSMC6wBpjyow6deqgqsyYMYOvvvrK6XCMyaNMJXugLvCdiPyIawKe/6jqZ8B0oLeI7Ad6u98bY0yZcf78eRYvXszQoUM5cOCA0+EY46VMN+MXhTXjG2NK24EDB2jfvj0hISFs2LCBoKAgp0MyFYxPNOMbY4wvu+WWW/joo4/Yt28f0dHRZGVlOR2SMYAle2OMKVZ33XUXr7zyCqtWrWLDhg1Oh2MMUPix8X3axYsXSUxMJD093elQjMlXQEAAoaGhNva6j3v00Ufp0aMHzZo1Izk5mSNHjhAREUHNmjWdDs1UUBUq2ScmJlK9enUiIiIQEafDMcaLqpKSkkJiYiINGzZ0OhxTBCJCs2bNAFi+fDljxowBICgoiPDwcMLDw5k3bx7169fn4MGDHDt2jPDwcG688Ub732RKRIVK9unp6ZboTZklItSpUwcbEKp86devHx9++CEJCQmer/j4eK677joA5s+fz7Rp0wBXy06DBg0IDw9nxYoVVKtWjR07dpCamkp4eDj169fHz8/PyY9jfFSFSvaAJXpTptnvZ/kTFhZGWFjYZcsfeeQROnbs6HUhcPToUc/FwEsvvcSSJUsA8PPzIywsjCZNmrBqlWui0W+//ZaMjAwiIiIICwujatWqJf+hjM+pcMneFyxatIjNmzfz6quvMm/ePKpVq8aIESMuu/6kSZPo1q0bvXr1uuw6sbGx+Pv707lz5yLFFh8fT//+/dm1a9cV13vhhRd45plnPO87d+7M+vXri3Tssqhfv368++67hb4Xm/NnbCq27Gb9y5k6dSrR0dHEx8d7LghyPjI9efJk1q5d63kfEhLCnXfeyQcffABATEwMlStX9hwnMDCw5D6MKbMs2Zdx2ff6rmTq1KlXXSc2NpbAwMAiJ/uCyp3sy1uiV1VU1VO7MqakXO1iYPHixRw4cMDrNkHOoXufeuop9u/f73lfu3Zthg4dyrx58wDXhWeNGjUIDw8nIiKCWrVqWQtTOWSP3pWyQYMG0bZtW5o1a8brr7/uWb5w4UJuu+02unfvzrp16zzLp0yZwsyZMwHYvn07nTp1IjIyksGDB3Py5EkARo0axUcffQRAREQEkydPpk2bNrRo0YK4uDji4+OZN28es2fPplWrVnz77bdeMX399de0atWKVq1a0bp1a06fPo2q8re//Y3mzZvTokULTy0hp0WLFjF27FjP+/79+xMbG8uECRM4f/48rVq1Ijo6GsBTm7jcfmNjY+nRowdDhgyhcePGREdHc7UBn0aNGsW4cePo3LkzN998s+ccxMbG0r9/f896Y8eOZdGiRZ7z88wzz3DHHXfQrl07tm7dSp8+fbjllls8//wA/ud//of27dsTGRnJ5MmTAVerRpMmTXj00Udp06YNhw8fJiIiwjO96ZIlS4iMjKRly5YMHz4cgJUrV9KxY0dat25Nr169+O233674mYwprLCwMHr06MHIkSOZNGkS8+fPZ8aMGZ7yr7/+mvXr1/Pee+/x4osvMmzYMBo3bgy4/h7//Oc/c88999C2bVvq1KlDUFCQ50JdVZk5cybvv/8+GzduJCkpiUuXLjnyOU3RVOiafY8ePfIsGzZsGI8++ijnzp2jX79+ecpHjRrFqFGjOH78OEOGDPEqi42NveoxFyxYQO3atTl//jzt27fn3nvv5cKFC0yePJktW7ZQo0YNevbsSevWrfNsO2LECObMmUP37t2ZNGkSzz33HC+//HKe9YKDg9m6dStz585l5syZvPnmm4wZM4bAwECeeuqpPOvPnDmTf//733Tp0oUzZ84QEBDAsmXL2L59Oz/++CPHjx+nffv2dOtWsAkIp0+fzquvvsr27dvzlF1pv9u2bWP37t3Uq1ePLl26sG7dOrp27XrFYyUlJfHdd98RFxdHVFRUnp9JfsLCwtiwYQNPPvkko0aNYt26daSnp9OsWTPGjBnD6tWr2b9/P5s2bUJViYqK4ptvvqFBgwbs27ePhQsXMnfuXK997t69m2nTprFu3TqCg4M5ceIEAF27dmXjxo2ICG+++SYzZszgpZdeKtB5NKY4hISEEBISwh133JFveWJiotctgoSEBFq0aAHAiRMn+Nvf/ua1vr+/Py+++CLjx48nLS2NV155xdP6EB4eTmhoKJUrV+jUUibZT6SU/etf/2L58uUAHD58mP3793P06FF69OjBDTfcAMB9993HTz/95LVdamoqp06donv37gCMHDmSoUOH5nuMe+65B4C2bduybNmyq8bUpUsXxo8fT3R0NPfccw+hoaF89913/PGPf8TPz4+6devSvXt3fvjhByIjI6/5swOX3W9QUBAdOnQgNDQUgFatWhEfH3/VZD9o0CAqVapE06ZNC1xrjoqKAqBFixacOXOG6tWrU716dQICAjh16hSrV69m9erVnguuM2fOsH//fk8v6U6dOuXZ51dffcWQIUMIDg4GXE2l4PpHet9995GUlMSFCxfskTpTpmQ/AVKnTh3atm2bp7xOnTqkpaXleZKgVatWgKu1a9KkSV7b+Pn5sXDhQoYPH05CQgJvvPGG50IgIiKCBg0aEBAQUBofz+RQoZP9lWri1apVu2J5cHBwgWryuY/35ZdfsmHDBqpVq0aPHj08A/wU5z2y7N64fn5+ZGZmXnX9CRMm8Ic//IFVq1bRqVMnvvzyy6s2oQNUrlzZq0mvIIMVXWm/OXsRFzT2nNtk7/tqcWVvU6lSJa/tK1WqRGZmJqrKxIkT+a//+i+v7eLj47n++uvzjUNV8/0ZPvbYY4wfP56oqChiY2OZMmXKVT+TMWVJ9erVad68Oc2bN89TFhkZyfnz5zl8+LBX60DLli0B2LdvH9OnT88zbPDq1avp3bs3W7Zs4f333yciIsKrdcDmFCh+ds++FKWmplKrVi2qVatGXFwcGzduBKBjx47ExsaSkpLCxYsX+fDDD/NsW6NGDWrVquW53/7WW295avkFUb16dU6fPp1v2YEDB2jRogVPP/007dq1Iy4ujm7duvHBBx+QlZXFsWPH+Oabb+jQoYPXdhEREWzfvp1Lly5x+PBhNm3a5CmrUqUKFy9ezHOsguw3t4kTJ3paQwoiPDycPXv2kJGRQWpqKmvWrCnwtgB9+vRhwYIFnDlzBoBff/2V5OTkK27zu9/9jqVLl5KSkgLgacZPTU2lfv36gKsjlTHlTUBAAI0aNaJ37948/PDDPP/8854WwLvvvpv09HTi4+P5+uuvWbJkCc8//zxNmjQBYM+ePcyZM4exY8cyYMAAIiMjqVGjhudpn88++4zHH3+cWbNm8fHHH7N582aOHz9eoMqI8Vaha/alrW/fvsybN4/IyEhuv/12T3NwSEgIU6ZM4Y477iAkJIQ2bdp4XQln1xgXL17MmDFjOHfuHDfffDMLFy4s8LEHDBjAkCFD+OSTT5gzZw533nmnp+zll19m7dq1+Pn50bRpU37/+9/j7+/Phg0baNmyJSLCjBkzuOmmm4iPj/ds16VLFxo2bEiLFi1o3rw5bdq08ZSNHj2ayMhI2rRpwzvvvONZPnjw4Hz3GxcXd9nYd+7c6Wl6L4iwsDCGDRtGZGQkjRo1yrf/w5Xcfffd7N2713OPMzAwkLfffvuKg5k0a9aMZ599lu7du+Pn50fr1q1ZtGgRU6ZMYejQodSvX59OnTpx8ODBQsVijK/L+dhfbsOHDyc6Oprk5GTPLYKEhAQiIiIA18VAzgvvbCdPnqRmzZosXryYDRs2eG4RZB+nXr169kRBLhVqitu9e/d6rih9xWOPPUabNm144IEHnA7FMX369OHzzz93OoxS44u/p8aUFFXl5MmTnlsEiYmJnqeAnn32WebNm+dpSQPXxXlaWhoiwvTp0/npp5+8bhE0bNjQczFRHl1uilur2Zdh//jHP/j+++8r/H3eipTojTHeRITatWtTu3btPK1006ZNY9q0aZw5c8ZzMZCamuqp1SckJPDZZ5+RlJTk2aZJkybs2bMHgHHjxpGSkuLVgbBRo0bcfPPNpfcBS4nP1OxFpC/wCuAHvKmq06+0fnmp2ZuKx35PjSleGRkZHD58mISEBDIzM+nTpw8A0dHRrF+/nsOHD3tund59992eCsbgwYPx8/PzukXQrFkzbrnlFsc+y9X4dM1eRPyAfwO9gUTgBxGJUdU9zkZmjDGmrKtatSq33nort956q9fy7P5EWVlZHDlyhISEBK/ppTMzM9mzZw//+c9/PE/1PPTQQ7z55ptcunSJDh06UK9ePa/bBG3bti2TLQM+keyBDsDPqvoLgIi8DwwECp3sL/eIlDFlga+0tBlTnmRPMJR7wqKVK1cCrr/L7E6E1atXB+DcuXOEhIR4njRIS0sDXKOeTp48meTkZLp16+Z1iyA8PJyuXbtecfjjkuIryb4+cDjH+0SgY2F3EhAQQEpKCnXq1LGEb8qc7PnsbcARY8oWEaFu3brUrVvXsywwMNBzMQBw6tQpEhISPANqZWRk0Lx5cxISEti2bZtn6ur58+fz4IMPsm3bNtLS0gr1CHVR+Eqyzy8z56kCichoYDRAgwYN8mwQGhpKYmKizRduyqyAgADPKILGGN9Rs2ZNr9kvw8LCPPN1gKsl4NChQ56RUv39/Uv13r+vJPtEIGf7SihwJPdKqvo68Dq4OujlLq9SpYoNV2qMMabUVatWzTMBEbjG5ihNvjKC3g9AIxFpKCL+wP1AjMMxGWOMMT7BJ2r2qpopImOBz3E9erdAVXc7HJYxxhjjE3wi2QOo6ipgldNxGGOMMb7GZwbVKSwROQYkFOMug4Hjxbi/isrOY9HZOSw6O4dFZ+ew6EriHIar6g25F5bbZF/cRGRzfqMSmcKx81h0dg6Lzs5h0dk5LLrSPIe+0kHPGGOMMdfIkr0xxhhTzlmyL7jXnQ6gnLDzWHR2DovOzmHR2TksulI7h3bP3hhjjCnnrGZvjDHGlHOW7AtARPqKyD4R+VlEJjgdj68RkTARWSsie0Vkt4g87nRMvkpE/ERkm4h86nQsvkhEaorIRyIS5/59vMPpmHyRiDzp/lveJSLviYjN3nQVIrJARJJFZFeOZbVF5AsR2e/+Xqukjm/J/ipExA/4N/B7oCnwRxFp6mxUPicT+KuqNgE6AX+xc3jNHgf2Oh2ED3sF+ExVGwMtsXNZaCJSHxgHtFPV5rhGNb3f2ah8wiKgb65lE4A1qtoIWON+XyIs2V9dB+BnVf1FVS8A7wMDHY7Jp6hqkqpudb8+jesfbH1no/I9IhIK/AF40+lYfJGIBAHdgPkAqnpBVU85GpTvqgxcJyKVgWrkMzGZ8aaq3wAnci0eCCx2v14MDCqp41uyv7r6wOEc7xOxRHXNRCQCaA1873Aovuhl4O/AJYfj8FU3A8eAhe5bIW+KyPVOB+VrVPVXYCZwCEgCUlV1tbNR+ay6qpoErkoRcGNJHciS/dVJPsvsEYZrICKBwMfAE6qa5nQ8vkRE+gPJqrrF6Vh8WGWgDfCaqrYGzlKCzabllfu+8kCgIVAPuF5E/uRsVOZqLNlfXSIQluN9KNZkVWgiUgVXon9HVZc5HY8P6gJEiUg8rltJd4nI286G5HMSgURVzW5V+ghX8jeF0ws4qKrHVPUisAzo7HBMvuo3EQkBcH9PLqkDWbK/uh+ARiLSUET8cXVEiXE4Jp8iIoLrPuleVZ3ldDy+SFUnqmqoqkbg+h38SlWtNlUIqnoUOCwit7sX/Q7Y42BIvuoQ0ElEqrn/tn+HdXS8VjHASPfrkcAnJXUgn5ni1imqmikiY4HPcfU6XaCqux0Oy9d0AYYDO0Vku3vZM+5pi40pTY8B77gv3H8BHnA4Hp+jqt+LyEfAVlxP2mzDRtO7KhF5D+gBBItIIjAZmA4sFZGHcF1EDS2x49sIesYYY0z5Zs34xhhjTDlnyd4YY4wp5yzZG2OMMeWcJXtjjDGmnLNkb4wxxpRzluyNMcaYcs6SvTHGGFPOWbI3xlwTEblORL4WET8Ricg5T3ch9+MvIt+4Z1AzxpQAS/bGmGv1ILBMVbOKshP31NFrgPuKJSpjTB6W7I0xeYjIWhHp7X79TxH5Vz6rRZPPWN4icrN7Ctn27hp/nHs62V0i8o6I9BKRdSKyX0Q6uDdb4d6fMaYEWLI3xuRnMvCsiEQDrYEncxa6x5a/WVXjcy2/Hdfshg+o6g/uxbcCrwCRQGPg/wBdgaeAZ9zr7ALal8gnMcbYRDjGmLxU9Rv3jGbjgR75NNUHA6dyLbsBV03/3lyTRR1U1Z0AIrIbWKOqKiI7gQj38bJE5IKIVFfV08X/iYyp2Kxmb4zJQ0RaACFAxmWS73kgINeyVOAwrlkOc8rI8fpSjveX8K5wVAXSrzVmY8zlWbI3xngRkRDgHWAgcFZE+uReR1VPAn4ikjPhXwAGASNE5P8U8ph1gGOqevGaAzfGXJYle2OMh4hUA5YBf1XVvcDzwJTLrL4a1713D1U9C/QHnhSRgYU4dE9gVaEDNsYUiM1nb4y5JiLSGhivqsOLYV/LgImquq/okRljcrOavTHmmqjqNmCtiPgVZT/unv0rLNEbU3KsZm+MMcaUc1azN8YYY8o5S/bGGGNMOWfJ3hhjjCnnLNkbY4wx5Zwle2OMMaacs2RvjDHGlHOW7I0xxphy7v8Bv9nafI9thrwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "M0 = flopy.modflow.Modflow(modelname='model', exe_name='../mf2005.exe')\n",
    "flopy.modflow.ModflowDis(M0, nlay=1, nrow=1, ncol=ncol, nper=1, delr=1., delc=1., top=0., botm=-b, steady=True, \n",
    "                         perlen=1., nstp=1)\n",
    "flopy.modflow.ModflowBas(M0, ibound=np.hstack([np.ones([nrow, ncol-1], dtype=int), -1*np.ones([1,1])]), \n",
    "                         strt=np.zeros([nrow, ncol], dtype=float))\n",
    "flopy.modflow.ModflowLpf(M0, hk=K, vka=-999., ss=-999., sy=-999., ipakcb=53)\n",
    "flopy.modflow.ModflowWel(M0, stress_period_data={0: [0, 0, ocol, 1.]})\n",
    "flopy.modflow.ModflowPcg(M0, hclose=1e-6, rclose=1e-6)\n",
    "flopy.modflow.ModflowOc(M0, stress_period_data={(0,0): ['save head', 'save budget']})\n",
    "M0.write_input()\n",
    "success, buff = M0.run_model(silent=True)\n",
    "A = flopy.utils.binaryfile.HeadFile('model.hds').get_data()[0,0,:]\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "plt.subplots(figsize=[8,2])\n",
    "plt.plot(X/1000., A, 'k--', mfc='none', label='adjoint solution, numerical')\n",
    "plt.xlabel('$x$ (km)')\n",
    "plt.ylabel(r'$\\psi^*_1$ (T/$L^2$)')\n",
    "plt.legend(loc=3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgAAAACdCAYAAADCH8BOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnjUlEQVR4nO3deVxV1fr48c8jWGqAY2aKiCmmEoribIZThUMOOYemDQ5lptbV9Jbmq9tgwy/TrkPaTcrsh9nNwpxKu06V3vLmkHNeUUlToVBwBHm+f5zjuaCAoMABzvN+vXhxzlp7n/2c7ZH1nLXWXltUFWOMMcZ4lhLuDsAYY4wxBc8SAGOMMcYDWQJgjDHGeCBLAIwxxhgPZAmAMcYY44EsATDGGGM8kLe7AyhIlSpV0sDAQHeHYYwxxhSILVu2xKvqrZnVFekEQEQigOmAF/C+qk7NbvvAwEB++umnAonNGGOMcTcROZRVXZEdAhARL2Am0AmoDwwQkfrujcoYY4wpGopsAgA0A35V1f+q6kUgGuju5piMMcaYIqEoDwFUA46kex4HNHdTLMYYU2xdunSJlJQU0tLSMvyUK1cOgNOnT3Pu3LkMdQDVq1cH4NixYyQlJWWo9/b2pm7dugDs37+fxMTEDPWlSpUiLCwMgJ9++umq+rJly9K6dWsA1qxZk6FeValcuTLt27cHYPHixZw+fTrD/jVr1iQiIgKAWbNmcfbs2Qz1d911F926dQNgypQpV73/Vq1a0bNnT1JTUxkzZsxV56ZTp0706tWL5ORkhg8fflX9Qw89RK9evTh+/DiPPPIIaWlp3HHHHcyaNatg/lEp2gmAZFJ21Y0NRGQYMAwgICDgqh1SUlKIi4vj/PnzeR6gMQWtVKlS+Pv7U7JkSXeHYoqRAQMGsHjx4gxl3t7epKSkADB69GiioqIy1JcvX54//vgDgKeffprPPvssQ31AQACHDjmGp0eNGsWqVasy1NevX5+dO3e6Xv/777/PUN+8eXM2bdoEwNixY9mxY0eG+o4dO7oSgPHjxxMbG5uhvmfPnq4EYPLkySQkJGSof/jhh10JwBtvvEFKSgolSpRw/Vy6dImePXsCEB0dnaGuRIkS1KpVC3AkT//+97+vqr98PFUlPj6eEiVKkJycTEGSonozIBFpCUxR1fudzycCqOprWe3TpEkTvXIS4MGDB/H19aVixYqIZJZTGFM0qCoJCQkkJSVRs2ZNd4djiqhz587x6aefMnv2bBYsWEBQUBDbtm1jxYoVGRowLy8vRo8eDcC3337L3r17M9TffPPNDBw4EIANGzZw5MiRDPW33HILnTp1AmDz5s2uRvDyj4+PDy1btgRg69atJCcnX1Vfv75j2tfevXu5ePEiJUqUQERcr3+5B+K3334jLS0tw/6lSpWibNmygKMHQ0Tw8vLK8P68vLwK9NznBxHZoqpNMq0rwgmAN7AP6AD8BvwIPKSqO7PaJ7MEYPfu3dStW9caf1MsqCp79uyhXr167g7FFDEHDhxgzpw5fPDBB/zxxx/Uq1eP9957jzZt2rg7NHMDsksAiuwQgKqmishTwCoclwF+kF3jnx1r/E1xYZ9lcz2SkpIICQnh4sWL9OzZk5EjRxIeHm6fp2KuKF8FgKouV9U6qlpLVV9xdzz5KSoqiqeeegqAOXPm8NFHH2W7/eTJk1m9enW226xdu/aqcbWCEBMTw9SpjiUbvvjiC3bt2uWqy0ncnTt3JjExkcTExAKdMJMf0p+L3AoMDCQ+Pj6PIzKe4OTJk0ydOpV+/foB4Ovry8KFCzl8+DCLFy+mbdu21vh7gCLbA+DJRowYcc1tXnrppWtus3btWnx8fGjVqlVehJVj3bp1c02u+eKLL+jatatrLC8ncS9fvhyA2NhYZs2axZNPPpl/weaj1NTUDOfCmPykqmzatImZM2eyePFiLl68SPv27UlOTsbHx8c1oc14jiLdA1Bc9OjRg7CwMIKDg5k7d66rfP78+dSpU4fw8HC+++47V/mUKVN46623AMfkmBYtWtCgQQN69uzJn3/+CcCQIUNcs24DAwN58cUXady4MSEhIezZs4fY2FjmzJnDtGnTCA0NZcOGDRliWrduHaGhoYSGhtKoUSOSkpIAePPNN2natCkNGjTgxRdfBBwNcb169Rg6dCjBwcHcd999nDt3DoAZM2ZQv359GjRoQP/+/YH/9WZ8//33xMTEMG7cOEJDQzlw4IAr7hUrVtC3b19XPGvXruWBBx5wvZ/4+HgmTJjAgQMHCA0NZdy4cQwaNIgvv/zStU9kZCQxMTFZnvfs4m7btq1r1cj4+HguLyEdFRVFjx49eOCBB6hZsyZ///vfefvtt2nUqBEtWrRwzXo+cOAAERERhIWF0aZNG/bs2eP6d3nmmWdo164dzz33XIaenePHj9OzZ08aNmxIw4YNXb0zWX0+jMmNRYsW0apVK5YuXcrw4cPZtWsXa9aswcfHx92hGXdRVY/5CQsL0yvt2rXrqrKClpCQoKqqZ8+e1eDgYI2Pj9ejR49q9erV9cSJE3rhwgVt1aqVjhw5UlVVX3zxRX3zzTdVVTUkJETXrl2rqqqTJk3S0aNHq6rq4MGDdfHixaqqWqNGDZ0xY4aqqs6cOVMfe+yxq17nSl27dtWNGzeqqmpSUpKmpKToqlWrdOjQoZqWlqaXLl3SLl266Lp16/TgwYPq5eWlP//8s6qq9unTRxcsWKCqqrfffrueP39eVVX//PNPVVWdP3++672kjzP985SUFK1evbomJyerquqIESNcr1mjRg09efKkHjx4UIODg137rl27Vrt3766qqomJiRoYGKgpKSlZnvfs4g4PD9cff/xRVVVPnjypNWrUcMVeq1YtPX36tJ44cUL9/Px09uzZqqo6ZswYnTZtmqqqtm/fXvft26eqqps2bdJ27dq53l+XLl00NTX1qnPRt29f1/6pqamamJioqpl/PtKfhysVhs+0cb89e/bo6NGjNSoqSlUd/4/fe+89TUpKcnNkpiABP2kWbaINAVyhbdu2V5X17duXJ598krNnz9K5c+er6ocMGcKQIUOIj4+nd+/eGerWrl17zWPOmDGDJUuWAHDkyBH279/P77//Ttu2bbn1Vsc9HPr168e+ffsy7Hfq1CkSExMJDw8HYPDgwfTp0yfTYzz44IMAhIWF8fnnn18zptatW/PMM88QGRnJgw8+iL+/P19//TVff/01jRo1AiA5OZn9+/cTEBBAzZo1CQ0NdR3j8jW3DRo0IDIykh49etCjR49rHvcyb29vIiIiWLp0Kb1792bZsmW88cYb2e4THh7OyJEjOXHiBJ9//jm9evXC2zv7j3hWcWenXbt2+Pr64uvrS9myZV09EyEhIWzfvp3k5GS+//77DP8WFy5ccD3u06dPppcXffvtt665HV5eXq5LlDL7fFSsWPGacRrPk5qaSkxMDLNmzWLNmjWULFmS8ePHA+Dj48OwYcPcHKEpTCwBcLO1a9eyevVqfvjhB8qUKUPbtm1dixLl5SScm2++GXA0LKmpqdfcfsKECXTp0oXly5fTokULVq9ejaoyceJEhg8fnmHb2NhY1+tfPsblrvRly5axfv16YmJi+Nvf/uZa2CMn+vXrx8yZM6lQoQJNmzbF19f3mvsMGjSIhQsXEh0dzQcffHDN7bOK29vb27Wa2ZWLRKXf5/L1zpcfp6amulZI27p1a6bHvOWWW64Z12XZfT6MuVLfvn1ZsmQJAQEBvPrqqzz22GNUrlzZ3WGZQsoSgCtk9429TJky2dZXqlQpR9/40zt16hTly5enTJky7Nmzx7WyVfPmzRk9ejQJCQn4+fmxePFiGjZsmGHfsmXLUr58eTZs2ECbNm1YsGCBqzcgJ3x9fTl9+nSmdQcOHCAkJISQkBB++OEH9uzZw/3338+kSZOIjIzEx8eH3377LdsV59LS0jhy5Ajt2rXj7rvv5pNPPrlqpStfX1/X/IIrtW3blscee4x58+a5Zitfa98hQ4bQrFkzqlSpQnBwMOBYBOThhx9mzZo12Z6P9AIDA9myZQvNmjW7agWza/Hz86NmzZosXryYPn36oKps3779qn+/K3Xo0IHZs2czZswYLl26xJkzZ7L8fBijqmzYsIE5c+Ywffp0br31VkaNGsWQIUPo0qVLsVjExuQvmwToZhEREaSmptKgQQMmTZpEixYtALj99tuZMmUKLVu2pGPHjjRu3DjDfpd7Bz788EPGjRtHgwYN2Lp1K5MnT87xsR944AGWLFmS6STAd955h7vuuouGDRtSunRpOnXqxH333cdDDz1Ey5YtCQkJoXfv3lk23uBYAnPgwIGEhITQqFEjxo4d61o7/LL+/fvz5ptv0qhRIw4cOJChzsvLi65du7JixQq6du161etXrFiR1q1bc9dddzFu3DgAbrvtNurVq8cjjzzi2u7YsWPXHAq40l/+8hdmz55Nq1atrutSu4ULF/KPf/yDhg0bEhwcnGFyYlamT5/Ov/71L0JCQggLC2Pnzp1Zfj6M50pKSmLWrFmEhIQQHh7OihUr2L59O+AYnurWrZs1/iZHiuxKgNcjq5UAi9qqaaNGjaJx48YZGjnjcPbsWUJCQvjPf/7jGkP/+9//TkBAgMdcblcUP9MmZ/78808CAwM5ffo0jRs3ZuTIkfTv358yZcq4OzRTSBXLlQA91aRJk9i8eTNTpkxxdyiFzurVq3n00Ud55plnXI0/4LrMzpii5uLFiyxZsoSdO3fy0ksvUb58eSZPnkybNm1o2rSpLdZjboj1ANi3JVPM2Ge66IuLi2Pu3LnMmzeP33//ndq1a7Njxw5KlSrl7tBMEZNdD4DNATDGmEJk0aJFBAYG8vLLLxMWFsby5cvZu3evNf4mz9kQAI7ZtNaVZooDT+rRKy4SExP58MMPqV+/Pvfeey/33HMPzz77LCNGjLDbOpt85fE9AKVKlSIhIcH+cJoiT1VJSEiwb4pFxNatWxk2bBjVqlVjzJgxLF26FHBcAfT6669b42/yXaHuARCRKcBQ4KSz6K+qutxZNxF4DLgEPK2qq67nGP7+/sTFxXHy5Mlrb2xMIVeqVCn8/f3dHYa5hkGDBvHxxx9TunRpIiMjeeKJJ6661NeY/FaoEwCnaar6VvoCEakP9AeCgarAahGpo6qXcvviJUuWtEzbGJOvVJWvvvqK++67j5tvvpmIiAiaNGnC4MGDr1obw5iCUlSHALoD0ap6QVUPAr8CzdwckzHGXGX9+vW0atWKbt268fHHHwOOO1WOHj3aGn/jVkUhAXhKRLaLyAciUt5ZVg04km6bOGeZMcYUCjt27KBr166Eh4dz5MgR3n//fQYPHuzusIxxcfsQgIisBqpkUvU8MBv4G6DO3/8PeBTIbMp+prP4RGQYMAwgICAgDyI2xpjsqSpDhw5l7969TJ06lVGjRtlqfabQcXsCoKodc7KdiMwDvnI+jQOqp6v2B45m8fpzgbngWAjo+iM1xpisJSQk8MYbbzBu3DgqVapEVFQUlStXpkKFCu4OzZhMFeohABG5Pd3TnsAvzscxQH8RuVlEagJBwL8LOj5jjDlz5gyvvvoqd9xxB2+99RarV68GoG7dutb4m0LN7T0A1/CGiITi6N6PBYYDqOpOEfkU2AWkAiOv5woAY4y5XqrKvHnzmDJlCseOHaNbt268+uqrrttQG1PYFeoEQFUHZVP3CvBKAYZjjDEuIsLXX39NzZo1+fTTT7n77rvdHZIxuVKohwCMMaYwWbduHXfffTe7d+8GICoqio0bN1rjb4okSwCMMeYatm/fTpcuXWjbti2xsbH89ttvAPj4+Nh9REyRZQmAMcZkQVUZNmwYoaGhfP/997zxxhvs37+fjh1zdPGSMYVaoZ4DYIwx7nDq1CnKli2LiFC5cmXGjRvHhAkTKF++/LV3NqaIsB4AY4xxOnPmDK+88goBAQF88803ALz88su8/vrr1vibYscSAGOMx0tJSeG9996jdu3avPDCC7Rr185WDjXFng0BGGM8mqrSvn1712z+f/7zn7Rq1crdYRmT76wHwBjjkTZu3EhqaioiwlNPPUVMTIzrzn3GeAJLAIwxHmXbtm106tSJNm3asHDhQgD69evHAw88YJf0GY9yzQRARLwKIhBjjMlPsbGxDBo0iEaNGrF582befPNN+vbt6+6wjHGbnMwBeB94JL8DMcaY/KKq9O7dm507dzJ+/Hiee+45m9VvPF5OhgC8ROSlKwtFxFdEFudDTMYYc8POnDnD66+/zqlTpxAR5s6dy/79+5k6dao1/saQswTgUaCliDx+uUBEGgI/AbvzKzBjjLkeKSkpzJkzh9q1azNhwgRiYmIAaNy4Mf7+/m6OzpjC45oJgKqmAg8CT4hIZxEZCnyJ4xa8k280ABHpIyI7RSRNRJpcUTdRRH4Vkb0icn+68jAR2eGsmyE2c8cYj6eqLF68mODgYJ544glq167Nd999x6BBWd5U1BiPlpNJgLOASOANYD7QD2iuqqvzKIZfcCQY6684bn2gPxAMRACz0k1InA0MA4KcPxF5FIsxpihYuBACA6FECcfvhQsREebPn89NN93E0qVL7ZI+Y64hJ5MAtwINgBCgJFAHmCciO4Adqhp9IwGo6m4gs8tvugPRqnoBOCgivwLNRCQW8FPVH5z7fQT0AFbcSBzGmCJi4UIYNgzOnnU8P3SItKFDKQEsWLCAcuXK4eVlFy8Zcy3XTABUdW765yLiz/8Sgi7ADSUA2agGbEr3PM5ZluJ8fGW5McYD6F//ilxu/J1KnDsHzz9PxchIN0VlTNGT66WAVTUOR6O7PKf7iMhqoEomVc+r6pdZ7ZbZ4bMpz+rYw3AMF9ja3sYUcZcuXUIOH870jwCHDxd0OMYUaTlKAETkZmA80BI4A/wH+FRVD+Rkf1W9nptnxwHV0z33B446y/0zKc/q2HOBuQBNmjTJMlEwxhROx48fJyYmhqFDh+Ll5UVi2bKUO3Xq6g0twTcmV3K6FPCbQA1gIo4GuCIwX0QG5ldgQAzQX0RuFpGaOCb7/VtVjwFJItLCOfv/YRxXJRhjipHt27fz6KOPEhAQwLBhw9i7dy8A5WbOhDJlMm5cpgy88oobojSm6MppAtAQ+IuqbgOaqepfgHuBkTcagIj0FJE4HL0Ly0RkFYCq7gQ+BXYBK3FcdnjJudsTOFYo/BU4gE0ANKbYOHjwIB06dKBhw4YsWrSIxx9/nL1793LnnXc6NoiMhLlzoUYNEHH8njvXUW6MyTFRvXavuIh0Bb5X1T9E5Escje4eYLCqts7nGPNMkyZN9KeffnJ3GMaYKyQnJ3Po0CGCg4NJTk6mVatWREZGMnToUCpUqODu8IwpskRki6o2yazumnMARKQlsEz/lyn0AXrjmNTXJ8+iNMZ4nCNHjvDuu+8yb948qlSpws6dO/Hx8WHbtm12Zz5j8llOJgEOBmaKyD4cXfErVfWT/A3LGFOcbdu2jddee43PPvsMVaVXr16MHTuWEiUco5LW+BuT/3KyDsAIABGpC3QCokSkLPAvHAnBd+nG5o0xJlOpqamkpKRQunRpdu/ezYoVKxgzZgyjRo2iRo0a7g7PGI+T00mAqOoeVZ2mqhFAe2AjjiGAzfkVnDGm6EtMTOStt96iVq1aTJs2DYBevXoRFxfHW2+9ZY2/MW6S44WARKQ9jnsCJOJYv387jisDLuRPaMaYomz//v3MmDGD+fPnc+bMGdq2bUvTpk0BKFmyJCVLlnRzhMZ4ttysBPgxjsv+vHEsBdwDx416aud9WMaYom7MmDF88803DBgwgLFjxxIaGurukIwx6eQmAfhVVZc4Hy/Oj2CMMUXThQsXiI6O5t1332Xx4sXUrFmTadOm4evry+233+7u8IwxmcjJ7YA/EpExwA8i8mz+h2SMKSpOnDjBSy+9RI0aNRgyZAjnz5/n+PHjANSpU8caf2MKsZz0AHyIYyXAKsD9IvI0sM35s11VrTfAGA+UnJxM7dq1SUpKonPnzowZM4aOHTvaJXzGFBE5uQxwDbDm8nMR8Qbq40gKmmPDAcZ4hLS0NFauXMm6det4/fXX8fHxYcaMGbRo0YK6deu6OzxjTC7ldCngG7obYGFhSwEbk3tnzpzho48+Yvr06ezdu5dq1aqxfft2W6LXmCIgu6WAb/RugHb3DWOKse+++47q1avz5JNP4ufnxyeffMLBgwet8TemGMjpVQANge6qmigizVS1tbNXYC2wMN+iM8YUuB9//NF13X5ISAidOnXiySefpFWrVja+b0wxkpsegMvbxovI2zjuEZB2owGISB8R2SkiaSLSJF15oIicE5Gtzp856erCRGSHiPwqIjPE/ioZc0NSUlL47LPPaN26Nc2aNeOFF14AwM/Pj4ULF9K6dWtr/I0pZnKUAKjqV6r6h/NpH+AnwIe8uRvgL8CDwPpM6g6oaqjzZ0S68tnAMCDI+RORB3EY45FmzJhBlSpV6NOnD8eOHeOdd95h+fLl7g7LGJPPcrMQEACqehHIs7sBqupuyPndv0TkdsBPVX9wPv8Ix6qEK/IqJmOKq7S0NH744Qeio6N5/vnnqVKlChUrViQiIoIBAwbQqVMnvLy83B2mMaYA5DoBKGA1ReRn4DTwgqpuAKrhmIh4WZyzzBiTCVVly5YtREdHs2jRIuLi4ihVqhSdOnWic+fOREZGEhlp83mN8TQFkgCIyGocCwld6XlV/TKL3Y4BAaqaICJhwBciEgxk1lWQ5bWMIjIMx3ABAQEBuQvcmCJKVTlz5gw+Pj4cPXqUpk2bUrJkSe6//36mTp1Kt27d8PX1dXeYxhg3KpAEQFU7Xsc+F4ALzsdbROQAUAfHN37/dJv6A0ezeZ25wFxwrAOQ2ziMKUr27t3LokWLiI6OplatWixdupRq1arx5Zdf0qZNG8qXL+/uEI0xhUShHQIQkVuBP1T1kojcgWOy339V9Q8RSRKRFsBm4GHgXXfGaoy7ffjhh7zzzjts3boVEaFNmzZ0797dVd+tWzc3RmeMKYxyehlgvhGRniISh2OVwWUisspZdQ+wXUS2AZ8BI9JdifAE8D7wK3AAmwBoPMxvv/3G9OnTOXfuHIBrXH/atGkcOXKEdevW8fjjj7s5SmNMYZajpYCLC1sK2BRlJ06c4J///CfR0dFs2LABVWXlypXcf//9pKWlUaKE2/N5Y0whk91SwIV2CMAY45jMJyLs37+funXrkpaWRr169ZgyZQr9+vXjzjvvBLDG3xiTa5YAGFPIJCUlERMTQ3R0NAEBAcycOZPatWvz6quv0rlzZ+666y5blc8Yc8MsATCmkFi+fDkffPABy5Yt4/z58/j7+9OsWTPAsVDWc8895+YIjTHFiSUAxrjJhQsXWLNmDREREZQoUYLly5ezceNGHn/8cfr370/Lli2ta98Yk29sEqAxBSg1NZVvv/2W6OholixZQmJiIhs3bqR169acOnUKHx8fW4rXGJNnbBKgMYXAzp07adu2LfHx8fj5+dGjRw/69+/v6uYvW7asmyM0xngSSwCMyQeqyubNm4mOjqZGjRqMHTuWoKAgunbtSvfu3YmIiKBUqVLuDtMY48EsATAmDy1dupQVK1awYsUKYmNjuemmmxgxwnEn65tuuon58+e7OUJjjHGwGUbGXKdDhw6xYMECxo8f7yr76KOPWLBgAcHBwURFRXHixAmmT5/uxiiNMSZzNgnQmFzYuHEjc+fOZf369Rw6dAiAcuXKceDAASpUqEB8fDzlypXD29s614wx7pfdJEDrATAmE6rKrl27mD17NgMGDGDPnj0AHDx4kJUrV9KkSRNmzJjB1q1bSUhIoEKFCgBUqlTJGn9jTJFgf6mMSSc2NpZnn32W9evXEx8fD0DVqlU5fPgwdevWZcCAAQwcONBW4jPGFHmWABiPlJqaytatW1m3bh3r1q3j3nvvZdSoUZQrV46ff/6Zzp07Ex4eTnh4OHfccYerwbdv98aY4sLtf81E5E3gAeAijlv7PqKqic66icBjwCXgaVVd5SwPA6KA0sByYLR60mQGk2uXb6qjqjz44IOsWbOGpKQkAIKCgujQoQPgGM//73//685QjTGmQLg9AQC+ASaqaqqIvA5MBJ4TkfpAfyAYqAqsFpE6qnoJmA0MAzbhSAAigBVuid4USufOnWPz5s2sX7+edevWAbBmzRpEhAoVKjBw4EDuuece7rnnHqpWrermaI0xpuC5PQFQ1a/TPd0E9HY+7g5Eq+oF4KCI/Ao0E5FYwE9VfwAQkY+AHlgC4NHOnj1LmTJlAJg4cSJvv/02Fy9eRERo2LAhHTp0cPUC/OMf/3BztMYY435uTwCu8CiwyPm4Go6E4LI4Z1mK8/GV5caDnDp1io0bN7q+4W/ZsoW4uDhuu+02GjZsyNNPP014eDitW7emfPny7g7XGGMKnQJJAERkNVAlk6rnVfVL5zbPA6nAwsu7ZbK9ZlOe1bGH4RguICAgIBdRm8IkISGBkiVL4ufnx+eff06fPn1IS0ujZMmSNGvWjPHjx3N5Gkj//v3p37+/myM2xpjCrUASAFXtmF29iAwGugId0k3miwOqp9vMHzjqLPfPpDyrY88F5oJjIaBcB28K3NmzZ/nqq6/Yt28f+/bt4+eff+aXX35h3rx5PP744zRp0oRJkyYRHh5O8+bNXV3/xhhjcs7tKwGKSATwNhCuqifTlQcDnwDNcEwCXAMEqeolEfkRGAVsxjEJ8F1VXX6tY9lKgIWDqrJ7927279/Pvn37XL+7devGM888wx9//EHFihUB8Pf3Jzg4mHvuuYeePXtSr149N0dvjDFFR2G/HfDfgZuBb5zXWm9S1RGqulNEPgV24RgaGOm8AgDgCf53GeAKbAJgoZOamsqhQ4dcjfu+ffuoWbMmzz77LADNmzcnOTkZgFtvvZWgoCBKly4NQIUKFdi+fTu1atWyb/fGGJNP3J4AqGrtbOpeAV7JpPwn4K78jMtcW1paGkePHnU18qrquvNdo0aN+OWXX1zb+vn50bu34wIPESE6OprKlSsTFBREuXLlrnrtkJCQAnkPxhjjqdw+BFCQbAgg91SVhIQE9u3bx7Fjx+jVqxcAw4cP5+OPP+bs2bOubYOCgti3bx8AUVFRpKWlUadOHYKCgqhcubItn2uMMQWssA8BmELg9OnT7N+/n8aNGyMizJs3j/fff599+/aRmJgIOJbBPXfuHN7e3oSGhnLLLbcQFBTkauT9/f83N3PIkCHueSPGGGNyxBIAD3HhwgWSk5MpW7Ys3t7ebNy4kaioKNf4/PHjxwGIi4ujWrVqpKWl4efnx0MPPeRq5OvUqUOJEo4bSD7xxBPufDvGGGNukCUAhVRqaionT54kOTmZpKQk1+/Q0FCqVatGbGwsCxYscJVf3uaFF14gLCyMVatWMWzYMFd5SkoKAJs2baJ58+bExsaybNkygoKC6Nq1q6uRL1u2LODo4h8+fLg7T4Exxph8ZAlAHkhJSXE1tGXKlKFSpUqcP3+elStXZmi8k5OT6dChA+Hh4Rw9ejRDA325EX/ttdcYMmQI27Zto0mTq4dtPv74YyIjIzl8+DCTJ0+mVKlS+Pj44OPjg6+vr+sGN5UrV6Zdu3b4+vq66nx8fLjtttsAiIyMZODAgQV6nowxxhQelgBcp4sXL1K1alWSk5O5cOGCq3zChAm89tprnD17lp49e2bYR0QoU6YM4eHhlChRgt9//x0fHx+qVq3qasQDAwMBCAwMZNasWfj6+mZoxGvXdlw00bp1ay5evEjJkiUzja9Ro0ZERUVlGb9NyDPGGM9mVwFcJ1Xl6aef5pZbbsnwDTs0NJSwsDAuXbrE9u3bMzTepUuXdo2hG2OMMfktu6sALAEwxhhjiqnsEgD7OmqMMcZ4IEsAjDHGGA9kCYAxxhjjgTxqDoCInAQO5eFLVgLi8/D1TObsPBcMO88Fx851wbDzDDVU9dbMKjwqAchrIvJTVpMrTN6x81ww7DwXHDvXBcPOc/ZsCMAYY4zxQJYAGGOMMR7IEoAbM9fdAXgIO88Fw85zwbFzXTDsPGfD5gAYY4wxHsh6AIwxxhgPZAnAdRCRCBHZKyK/isgEd8dTXIlIdRH5l4jsFpGdIjLa3TEVZyLiJSI/i8hX7o6luBKRciLymYjscX6uW7o7puJIRMY6/2b8IiL/X0RKuTumwsgSgFwSES9gJtAJqA8MEJH67o2q2EoFnlXVekALYKSd63w1Gtjt7iCKuenASlWtCzTEzneeE5FqwNNAE1W9C/AC+rs3qsLJEoDcawb8qqr/VdWLQDTQ3c0xFUuqekxV/+N8nITjj2U190ZVPImIP9AFeN/dsRRXIuIH3AP8A0BVL6pqoluDKr68gdIi4g2UAY66OZ5CyRKA3KsGHEn3PA5rlPKdiAQCjYDNbg6luHoHGA+kuTmO4uwO4CQw3znU8r6I3OLuoIobVf0NeAs4DBwDTqnq1+6NqnCyBCD3JJMyu5QiH4mID/BPYIyqnnZ3PMWNiHQFTqjqFnfHUsx5A42B2araCDgD2ByiPCYi5XH0ytYEqgK3iMhA90ZVOFkCkHtxQPV0z/2x7qV8IyIlcTT+C1X1c3fHU0y1BrqJSCyOIa32IvKxe0MqluKAOFW93Iv1GY6EwOStjsBBVT2pqinA50ArN8dUKFkCkHs/AkEiUlNEbsIxuSTGzTEVSyIiOMZLd6vq2+6Op7hS1Ymq6q+qgTg+z9+qqn1jymOq+jtwRETudBZ1AHa5MaTi6jDQQkTKOP+GdMAmW2bK290BFDWqmioiTwGrcMwu/UBVd7o5rOKqNTAI2CEiW51lf1XV5e4LyZgbMgpY6Pzy8F/gETfHU+yo6mYR+Qz4D44riX7GVgTMlK0EaIwxxnggGwIwxhhjPJAlAMYYY4wHsgTAGGOM8UCWABhjjDEeyBIAY4wxxgNZAmCMMcZ4IEsAjDHGGA9kCYAxJs+ISGkRWSciXiISKCK/XOfr3CQi6513czPG5ANLAIwxeelR4HNVvXQjL+K81fYaoF+eRGWMuYolAMaYHBGRf4nIvc7HL4vIjEw2iwS+zGTfO5y3wG3q7BnY47wd7i8islBEOorIdyKyX0SaOXf7wvl6xph8YAmAMSanXgSeF5FIoBEwNn2lc337O1Q19oryO3Hc0fERVf3RWVwbmA40AOoCDwF3A38B/urc5hegab68E2OM3QzIGJMzqrreeXe1Z4C2mXTzVwISryi7FUePQK8rbpp1UFV3AIjITmCNqqqI7AACnce7JCIXRcRXVZPy/h0Z49msB8AYkyMiEgLcDlzIokE+B5S6ouwUcATHnR3Tu5DucVq652lk/GJyM3D+emM2xmTNEgBjzDWJyO3AQqA7cEZE7r9yG1X9E/ASkfRJwEWgB/CwiDyUy2NWBE6qasp1B26MyZIlAMaYbIlIGeBz4FlV3Q38DZiSxeZf4xjLd1HVM0BXYKyIdM/FodsBy3MdsDEmR0RV3R2DMaaYEJFGwDOqOigPXutzYKKq7r3xyIwxV7IeAGNMnlHVn4F/iYjXjbyO84qCL6zxNyb/WA+AMcYY44GsB8AYY4zxQJYAGGOMMR7IEgBjjDHGA1kCYIwxxnggSwCMMcYYD2QJgDHGGOOBLAEwxhhjPND/ARSrzHrmhQeLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#dhdK = np.empty(0)\n",
    "dhdK = []\n",
    "for oc,xp in enumerate(X[::1000]):\n",
    "    M1 = M0\n",
    "    flopy.modflow.ModflowWel(M1, stress_period_data={0: [0, 0, oc*1000, 1.]})\n",
    "    M1.write_input()\n",
    "    success, buff = M1.run_model(silent=True)\n",
    "    A = flopy.utils.binaryfile.HeadFile('model.hds').get_data()\n",
    "    #dhdK = np.append(dhdK, np.sum(A[:round(2*L/3)]*b*-R/K*b))\n",
    "    dhdK.append(-b*np.float(np.atleast_2d(np.gradient(np.ravel(H0)))@np.atleast_2d(np.gradient(np.ravel(A))).T))\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "f,s = plt.subplots(figsize=[8,2])\n",
    "plt.plot(X[::1000]/1000., dhdK, 'k--', mfc='none', label='adjoint sensitivity, numerical')\n",
    "plt.plot(ocol/1000., dhdK[int(ocol/1000)], 'ro')\n",
    "plt.xlabel('$x$ (km)')\n",
    "plt.ylabel(r'$\\partial h / \\partial K$')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src onerror=\"var cell = this.closest('.code_cell');var editor = cell.querySelector('.input_area');editor.style.background='rgba(0, 200, 0, 0.2)';this.parentNode.removeChild(this)\">"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-99.976749 \n"
     ]
    }
   ],
   "source": [
    "set_background('rgba(0, 200, 0, 0.2)')\n",
    "#print('%.6f (%.6f)'% (dhdK[int(np.where(X[::1000]==float(ocol))[0])], benchmark))\n",
    "print('%.6f '% (dhdK[int(np.where(X[::1000]==float(ocol))[0])]))"
   ]
  },
  {
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}