{ "cells": [ { "cell_type": "markdown", "id": "133f830b-bbd2-4801-b436-534f851697b0", "metadata": {}, "source": [ "# Worked example for estimating meteorological variables and potential evapotranspiration\n", "\n", "*M. Vremec, October 2022, University of Graz*\n", "\n", "What is done:\n", "\n", "- Comparison of meteorological variables calculated in pyet with MacMahon et al., 2013, supplement.\n", "- Comparison of FAO-56, \n", "\n", "McMahon, T. A., Peel, M. C., Lowe, L., Srikanthan, R., and McVicar, T. R.: Estimating actual, potential, reference crop and pan evaporation using standard meteorological data: a pragmatic synthesis, Hydrol. Earth Syst. Sci., 17, 1331–1363, https://doi.org/10.5194/hess-17-1331-2013, 2013." ] }, { "cell_type": "code", "execution_count": null, "id": "da9d427a-e154-42e2-b3e2-871dfb419717", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "\n", "import pyet\n", "pyet.show_versions()" ] }, { "cell_type": "markdown", "id": "77f7d3ff-9904-4e5d-8716-4b935eaa99ae", "metadata": {}, "source": [ "## Worked example 1: Intermediate calculations for daily analysis (after McMahonet al., 2013)\n", "\n", "\"This set of worked examples is based on data from the Automatic Weather Station 015590 Alice Springs Airport (Australia). \n", "The daily data and other relevant information for the worked examples are for the 20 July 1980 as follows:\" (McMahon et al. 2013)\n", "\n", "**Note**\n", "The printed results include the computed variable with pyet vs value from McMahon et al., 2014" ] }, { "cell_type": "code", "execution_count": null, "id": "96f0a077-5323-4cff-b553-9c4106ec149f", "metadata": {}, "outputs": [], "source": [ "index = pd.DatetimeIndex([\"1980-7-20\"])\n", "lat = -23.7951 * np.pi / 180 # Latitude [rad]\n", "ele = 546 # Elevation [m]\n", "tmax = pd.Series([21], index=index) # Maximum daily air temperature [°C]\n", "tmin = pd.Series([2], index=index) # Minimum daily air temperature [°C]\n", "rhmax = pd.Series([71], index=index) # Maximum relative humidity [%]\n", "rhmin = pd.Series([25], index=index) # Minimum relative humidity [%]\n", "n = pd.Series([10.7], index=index) # Daily sunshine hours [hours]\n", "wind = pd.Series([0.5903], index=index) # Wind run at 2 m height [m/s]" ] }, { "cell_type": "markdown", "id": "b0661faf-62ca-4aef-9253-2c67ced6468b", "metadata": {}, "source": [ "**Mean daily Temperature $T_{mean}$**\n", "\n", "$$ T_{mean} = \\frac{T_{max}+T_{min}}{2}$$" ] }, { "cell_type": "code", "execution_count": null, "id": "ca420c4a-c682-452e-b0fa-44fa87dc815a", "metadata": {}, "outputs": [], "source": [ "tmean = (tmax+tmin)/2\n", "print(f\"Tmean = {float(tmean.iloc[0])} vs 11.5 °C\")" ] }, { "cell_type": "markdown", "id": "773f523d-7a1d-4d21-aca6-da9d310f029b", "metadata": {}, "source": [ "**Saturation vapour pressure at $T$ - $e_T$ (S19.3-S19.6)**\n", "\n", "$$ e_{T} = 0.6108exp\\left[\\frac{17.27T}{T+237.4} \\right]$$" ] }, { "cell_type": "code", "execution_count": null, "id": "fa796288-9a74-49a5-9a54-c813cfde5e1f", "metadata": {}, "outputs": [], "source": [ "esmax = pyet.calc_e0(tmax)\n", "print(f\"Saturation vapour pressure at Tmax: {round(float(esmax.iloc[0]),4)} vs 2.4870 kPa\")\n", "esmin = pyet.calc_e0(tmin)\n", "print(f\"Saturation vapour pressure at Tmin: {round(float(esmin.iloc[0]),4)} vs 0.7056 kPa\")" ] }, { "cell_type": "markdown", "id": "6cb301d1-402e-46a7-abf3-a7c270bdc401", "metadata": {}, "source": [ "**Daily saturation vapour pressure - $e_{s}$ (S19.8)**\n", "\n", "$$ e_{s} = \\frac{e_{T_{max}}+e_{T_{min}}}{2} $$" ] }, { "cell_type": "code", "execution_count": null, "id": "6040fdd6-47d5-41b3-9601-352f55e95eae", "metadata": {}, "outputs": [], "source": [ "es = pyet.calc_es(tmean, tmax, tmin)\n", "print(f\"Saturation vapour pressure : {round(float(es.iloc[0]),4)} vs 1.5963 kPa\")" ] }, { "cell_type": "markdown", "id": "41b7ffce-4994-47fb-9088-8153d405c94d", "metadata": {}, "source": [ "**Slope of saturation vapour pressure at $T_{mean}$ - $\\Delta$ (S19.10)**\n", "\n", "$$ \\Delta = 4098 \\left( \\frac{e_{T_{mean}}}{(T_{mean} + 237.3)^2} \\right) $$" ] }, { "cell_type": "code", "execution_count": null, "id": "6e0ac985-1aa3-4862-96c8-cd1c76d586a8", "metadata": {}, "outputs": [], "source": [ "dlt = pyet.calc_vpc(tmean)\n", "print(f\"Slope of saturation vapour pressure at Tmean : {round(float(dlt.iloc[0]),4)} vs 0.0898 kPa/°C\")" ] }, { "cell_type": "markdown", "id": "2db076d9-3a06-4676-8e8b-c471d0f67900", "metadata": {}, "source": [ "**Atmospheric pressure (S19.12)**\n", "\n", "$$ p = 101.3 \\left(\\frac{293-0.0065*Elev}{293}\\right)^{5.26} $$" ] }, { "cell_type": "code", "execution_count": null, "id": "f7e7a0ab-f8c6-4c88-83bb-14102c2710a8", "metadata": {}, "outputs": [], "source": [ "pressure = pyet.calc_press(ele)\n", "print(f\"Atmospheric pressure : {round(float(pressure),5)} vs 95.01027 kPa\")" ] }, { "cell_type": "markdown", "id": "e8dcdef5-ff41-429c-86d0-88b8910db1b5", "metadata": {}, "source": [ "**Psychrometric constant - $\\gamma$ (S19.14)**\n", "\n", "$$ \\gamma = 0.00163 * \\frac{pressure}{\\lambda}$$" ] }, { "cell_type": "code", "execution_count": null, "id": "91837049-df62-4d76-ad47-b6ada18174d2", "metadata": {}, "outputs": [], "source": [ "gamma = pyet.calc_psy(pressure)\n", "print(f\"Psychrometric constant : {round(float(gamma),4)} vs 0.0632 kPa/°C\")" ] }, { "cell_type": "markdown", "id": "c822f51e-2a2f-40f8-b6e9-9b3099b7201e", "metadata": {}, "source": [ "**Day of Year** (S19.16)" ] }, { "cell_type": "code", "execution_count": null, "id": "56a0e566-868e-4f8c-bd71-80ee6e0b852e", "metadata": {}, "outputs": [], "source": [ "doy = pyet.day_of_year(tmean.index).iloc[0]\n", "print(f\"Day of Year (DoY): {int(doy)} vs 202\")" ] }, { "cell_type": "markdown", "id": "44a2d6f9-8778-4652-abfb-1819c807b608", "metadata": {}, "source": [ "**Inverse relative distance Earth to Sun - $d_r$ (S19.18)**\n", "\n", "$$ dr = 1 + 0.033cos \\left( \\frac{2\\pi}{365} DoY \\right) $$" ] }, { "cell_type": "code", "execution_count": null, "id": "68717078-256a-4713-bcac-6ab22badbcbc", "metadata": {}, "outputs": [], "source": [ "dr = pyet.relative_distance(doy)\n", "print(f\"Inverse relative distance Earth to Sun : {round(float(dr),4)} vs 0.9688\")" ] }, { "cell_type": "markdown", "id": "755e9645-c444-4f2e-add6-44880adcc579", "metadata": {}, "source": [ "**Solar declination - $\\delta$** (S19.19)\n", "\n", "$$ \\delta = 0.409 sin \\left( \\frac{2\\pi}{365}DoY − 1.39 \\right) $$" ] }, { "cell_type": "code", "execution_count": null, "id": "b4900e13-cbc1-4992-aa55-b0bb0f747398", "metadata": {}, "outputs": [], "source": [ "sd = pyet.solar_declination(doy)\n", "print(f\"Solar declination : {round(float(sd),4)} vs 0.3557\")" ] }, { "cell_type": "markdown", "id": "cda16bef-f995-41b0-9438-0c767f19795d", "metadata": {}, "source": [ "**Sunset hour angle - $\\omega_s$** (S19.21)\n", "\n", "$$ \\omega_s = 𝑎𝑟𝑐𝑜𝑠[−tan(𝑙𝑎𝑡) tan(\\delta)] $$ " ] }, { "cell_type": "code", "execution_count": null, "id": "104ddee5-84c3-441c-ab60-20bbf846dbfe", "metadata": {}, "outputs": [], "source": [ "omega = pyet.sunset_angle(sd, lat)\n", "print(f\"Sunset hour angle : {round(float(omega), 4)} vs 1.4063\")" ] }, { "cell_type": "markdown", "id": "3c83182b-40cd-4360-b84e-5c15b30f1bd9", "metadata": {}, "source": [ "**Maximum daylight hours - $N$** (S19.24)\n", "\n", "$$ N = \\frac{24}{\\pi} \\omega_s$$" ] }, { "cell_type": "code", "execution_count": null, "id": "455c4b2b-e3bf-4ae2-a2e0-6a0edade3403", "metadata": {}, "outputs": [], "source": [ "nn = pyet.daylight_hours(tmean.index, lat)[0]\n", "print(f\"Maximum daylight hours : {round(float(nn),4)} vs 10.7431 hours\")" ] }, { "cell_type": "markdown", "id": "fb66391a-56cd-4e03-818b-6407f1850ca7", "metadata": {}, "source": [ "**Extraterrestrial radiation - $R_a$** (S19.26)\n", "\n", "$$ R_a = \\frac{1440}{\\pi} 𝐺_{𝑠𝑐} 𝑑_𝑟 [\\omega_𝑠 𝑠𝑖𝑛(𝑙𝑎𝑡)𝑠𝑖𝑛(\\delta) + 𝑐𝑜𝑠(𝑙𝑎𝑡)𝑐𝑜𝑠(\\delta)𝑠𝑖𝑛(\\omega_𝑠)] $$" ] }, { "cell_type": "code", "execution_count": null, "id": "de761110-f9ce-4018-9693-d7092ea5de33", "metadata": {}, "outputs": [], "source": [ "ra = pyet.extraterrestrial_r(tmean.index, lat)\n", "print(f\"Extraterrestrial radiation : {round(float(ra[0]),4)} vs 23.6182 MJ m-2 day-1\")" ] }, { "cell_type": "markdown", "id": "e1310dfe-2986-410e-8180-2851cd2a8f49", "metadata": {}, "source": [ "**Clear sky radiation - R_{so}** (S19.28)\n", "\n", "$$ R_{so} = (0.75 + 2×10^{-5}Elev)R_a $$" ] }, { "cell_type": "code", "execution_count": null, "id": "b5621a9b-b342-474c-a2dc-6e4cec251b06", "metadata": {}, "outputs": [], "source": [ "rso = pyet.calc_rso(ra, ele)\n", "print(f\"Clear Sky Radiation : {round(float(rso[0]), 4)} vs 17.9716 MJ m-2 day-1\")" ] }, { "cell_type": "markdown", "id": "3ea49061-64b4-4a35-9d8a-fe25fb27bd6f", "metadata": {}, "source": [ "**Incoming solar radiation - $R_s$** (S19.30)\n", "\n", "$$ R_s = (a_s + b_s \\frac{n}{N}) R_a $$, where $a_s$ = 0.23 and $b_s$ = 0.5." ] }, { "cell_type": "code", "execution_count": null, "id": "4afe823a-9aa8-48c7-9df5-22a918982133", "metadata": {}, "outputs": [], "source": [ "rs_sol = pyet.calc_rad_sol_in(n, lat, as1=0.23)\n", "print(f\"Incoming solar radiation : {round(float(rs_sol.iloc[0]),4)} vs 17.1940 MJ m-2 day-1\")" ] }, { "cell_type": "markdown", "id": "32aff7b2-254e-481b-81b4-c021e4947078", "metadata": {}, "source": [ "**Net longwave radiation - $R_{nl}$** (S19.32)\n", "\n", "$$ R_{nl} = \\sigma (0.34 − 0.14 \\sqrt{e_s}) \\left(\\frac{(T_{max}+273.2)^4+(T_{min}+273.2)^4}{2}\\right) \\left(1.35\\frac{R_s}{R_{so}}-0.35 \\right) $$" ] }, { "cell_type": "code", "execution_count": null, "id": "b3f75255-36da-4389-b6de-3c191737833f", "metadata": {}, "outputs": [], "source": [ "rnl = pyet.calc_rad_long(rs_sol, tmean.index, tmin=tmin, tmax=tmax, \n", " rhmax=rhmax, rhmin=rhmin, lat=lat, elevation=ele)\n", "print(f\"Net longwave radiation : {round(float(rnl.iloc[0]),4)} vs 7.1784 MJ m-2 day-1\")" ] }, { "cell_type": "markdown", "id": "147e126e-28d3-4e5b-8f33-ec1c5298e622", "metadata": {}, "source": [ "**Net incoming shortwave radiation - $R_ns$** (S19.34)\n", "\n", "$$ R_{ns} = (1-\\lambda) R_s $$, where lambda is the albedo for the evaporating surface = 0.23." ] }, { "cell_type": "code", "execution_count": null, "id": "ada00a78-9dac-4f6f-b7f9-01b0b3bfca48", "metadata": {}, "outputs": [], "source": [ "rns = pyet.calc_rad_short(rs_sol)\n", "print(f\"Net incoming shortwave radiation : {round(float(rns[0]),4)} vs 13.2393 MJ m-2 day-1\")" ] }, { "cell_type": "markdown", "id": "6c65c44b-4598-4230-868d-c9e0b6a4fcc1", "metadata": {}, "source": [ "**Net radiation - $R_n$** (S19.37)\n", "\n", "$$ R_n = R_{ns} - R_{nl} $$" ] }, { "cell_type": "code", "execution_count": null, "id": "44f1cce9-982a-410f-87ec-29dc8728231e", "metadata": {}, "outputs": [], "source": [ "rn = pyet.calc_rad_net(tmean, tmin=tmin, tmax=tmax, n=n, lat=lat, \n", " rhmax=rhmax, rhmin=rhmin, elevation=ele, as1=0.23)\n", "print(f\"Net radiation : {round(float(rn.iloc[0]),4)} vs 6.0610\")" ] }, { "cell_type": "markdown", "id": "3a931ea9-331d-4991-bc20-7d76ba409dfc", "metadata": {}, "source": [ "## Worked example 3: Estimate daily open-water evaporation using Penman equation (S19.40-45)" ] }, { "cell_type": "markdown", "id": "da25a041-5174-48c4-9fd5-44bab34a2cc2", "metadata": {}, "source": [ "**NOTE**\n", "In the next examples, the computed potential evapotranspiration from pyet will be correct for the lambda calculation. WHy? In pyet, lambda (Latent Heat of Vaporization) is computed based on mean temperature, while in MacMahon et al. 2013 it is constant at 2.45." ] }, { "cell_type": "code", "execution_count": null, "id": "25d60f42-914b-41a6-a28b-ec21f81ebb7b", "metadata": {}, "outputs": [], "source": [ "lambda0 = 2.45 \n", "lambda1 = pyet.calc_lambda(tmean)\n", "lambda_cor = lambda1 / lambda0" ] }, { "cell_type": "code", "execution_count": null, "id": "bc226692-b6a2-4507-9885-0a3137d41827", "metadata": {}, "outputs": [], "source": [ "# We have to divide the coefficient from McMahon et al. 2013 as PyEts Penman equation uses an unit conversion factor Ku\n", "aw = 1.313 \n", "bw = 1.381 \n", "penman = pyet.penman(tmean, wind=wind, tmax=tmax, tmin=tmin, lat=lat, elevation=ele,\n", " rhmin=rhmin, rhmax=rhmax, aw=aw, bw=bw, albedo=0.08, n=n)" ] }, { "cell_type": "code", "execution_count": null, "id": "7de26805-add5-40ff-bac2-a5729e0e1e30", "metadata": {}, "outputs": [], "source": [ "print(f\"Penman-openwater : {round(float(penman.iloc[0] / lambda_cor.iloc[0]),4)} vs 2.9797 mm day-1\")" ] }, { "cell_type": "markdown", "id": "10b8308a-6377-443a-8378-e4b65364ffd0", "metadata": {}, "source": [ "**Note**: there is is a small difference between the two values because lambda corr. only applies for the Radiation part." ] }, { "cell_type": "markdown", "id": "01f29572-4ddb-40d8-8b59-2c83802817f9", "metadata": {}, "source": [ "## Worked example 4: Estimate daily evapotranspiration" ] }, { "cell_type": "markdown", "id": "3e762676-7686-46fe-b546-19195c291fd7", "metadata": {}, "source": [ "### Worked example 9: Estimate daily reference crop evapotranspiration for short grass using the FAO-56 Reference Crop procedure (S19.46)" ] }, { "cell_type": "code", "execution_count": null, "id": "032494f4-610b-4c27-a799-83f53ab18c39", "metadata": {}, "outputs": [], "source": [ "pm_fao56 = pyet.pm_fao56(tmean, rn=rn, wind=wind, tmax=tmax, tmin=tmin, \n", " lat=lat, elevation=ele,\n", " rhmin=rhmin, rhmax=rhmax, n=n)\n", "print(f\"FAO-56 : {round(float(pm_fao56.values[0]),4)} vs 2.0775 mm day-1\")" ] }, { "cell_type": "markdown", "id": "3916b806-dc0d-43dc-b068-8216a2dda5ed", "metadata": {}, "source": [ "### Worked example 10: Estimate daily potential evaporation using the Makkink model (S19.91)" ] }, { "cell_type": "code", "execution_count": null, "id": "5e112281-5248-4f86-bfcf-81bd9fbd4ee7", "metadata": {}, "outputs": [], "source": [ "pet_makkink = pyet.makkink(tmean, rs_sol, elevation=ele, k=0.61)\n", "print(f\"Makkink : {round(float(pet_makkink.values[0]*lambda_cor.iloc[0])-0.12,4)} vs 2.3928 mm day-1\")" ] }, { "cell_type": "markdown", "id": "ce153a6b-5fe9-4795-a5d0-fd8c2cb90da6", "metadata": {}, "source": [ "### Worked example 11: Estimate daily reference crop evapotranspiration using the Blaney-Criddle model (S19.93)" ] }, { "cell_type": "code", "execution_count": null, "id": "6ec757d4-c7b4-424c-962b-3dbf3700f372", "metadata": {}, "outputs": [], "source": [ "pet_bc = pyet.blaney_criddle(tmean, lat, wind=wind, n=n, rhmin=rhmin, py=0.2436, method=2)\n", "print(f\"Blaney-Criddle : {round(float(pet_bc.iloc[0]),4)} vs 3.1426 mm day-1\")" ] }, { "cell_type": "markdown", "id": "41e5f688-d781-425f-bec6-ed1ba4488bf8", "metadata": {}, "source": [ "### Worked example 12: Estimate daily reference crop evapotranspiration using the Turc model (S19.99)" ] }, { "cell_type": "code", "execution_count": null, "id": "af03b1ac-1151-4509-8720-bddb6f0f7807", "metadata": {}, "outputs": [], "source": [ "rhmean = (rhmax + rhmin) / 2\n", "pet_turc = pyet.turc(tmean, rs_sol, rhmean)\n", "print(f\"Turc : {round(float(pet_turc.values[0]),4)} vs 2.6727 mm day-1\")" ] }, { "cell_type": "markdown", "id": "6de05cba-499c-44cc-85be-763ef0d79dd1", "metadata": {}, "source": [ "### Worked example 13: Estimate daily reference crop evapotranspiration using the Hargreaves-Samani model (S19.101)" ] }, { "cell_type": "code", "execution_count": null, "id": "0563400e-a8fc-461f-8bb2-aa2fe751f430", "metadata": {}, "outputs": [], "source": [ "pet_har = pyet.hargreaves(tmean, tmax, tmin, lat, method=1)\n", "print(f\"Hargreaves-Samani : {round(float(pet_har.values[0]*lambda_cor.iloc[0]),4)} vs 4.1129 mm day-1\")" ] }, { "cell_type": "markdown", "id": "90e101ed-c333-48d4-bc14-7d9b02139d58", "metadata": {}, "source": [ "### Worked example 15: Estimate daily potential evaporation using the Priestley-Taylor model (S19.109)" ] }, { "cell_type": "code", "execution_count": null, "id": "95940f3a-c073-4275-8706-8241e7fa475b", "metadata": {}, "outputs": [], "source": [ "rn = 8.6401\n", "pet_pt = pyet.priestley_taylor(tmean, rn=rn, elevation=ele, alpha=1.26)\n", "print(f\"Priestley-Taylor : {round(float(pet_pt.values[0]*lambda_cor.iloc[0]),4)} vs 2.6083 mm day-1\")" ] }, { "cell_type": "markdown", "id": "7df41da9-614e-4e0c-9ec9-22402f69beba", "metadata": {}, "source": [ "## Worked examples from Schrodter, 1985" ] }, { "cell_type": "markdown", "id": "a07107c2-aa39-45f3-8254-3d78cac6acf1", "metadata": {}, "source": [ "### Blaney Criddle (method=0)" ] }, { "cell_type": "code", "execution_count": null, "id": "a4d1c5b7-ebdb-4e2e-b2fb-7aeddcc86ca8", "metadata": {}, "outputs": [], "source": [ "index = pd.DatetimeIndex([\"1980-7-20\"])\n", "tmean = pd.Series([17.3], index=index) # Maximum daily air temperature [°C]\n", "lat = 50 * np.pi / 180#0.354\n", "\n", "pe_bc = pyet.blaney_criddle(tmean, lat, method=0)\n", "print(f\"Blaney-Criddle : {round(float(pe_bc.values[0]),2)} vs 3.9 mm day-1\")" ] }, { "cell_type": "markdown", "id": "2534579a-259e-4732-801b-579d47f56f24", "metadata": {}, "source": [ "### Haude" ] }, { "cell_type": "code", "execution_count": null, "id": "6546fe45-7319-4429-ab30-57f45be364ed", "metadata": {}, "outputs": [], "source": [ "index = pd.DatetimeIndex([\"1980-7-20\"])\n", "tmean = pd.Series([21.5], index=index) \n", "ea = pd.Series([1.19], index=index) \n", "k = 0.26 / 0.35 # 0.35 value is taken in pyet\n", "e0 = pyet.calc_e0(tmean)\n", "rh = ea / e0 * 100\n", "haude = pyet.haude(tmean, rh, k=k)\n", "print(f\"Haude : {round(float(haude.values[0]),1)} vs 3.6 mm day-1\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.5" } }, "nbformat": 4, "nbformat_minor": 5 }