Hourly FAO-56 Penman-Monteith Reference Evapotranspiration (ET0) from EPW Files

Page 1 of 6 Method Statement Hourly FAO-56 Penman–Monteith Reference Evapotranspiration (𝑬𝑻𝟎) from EPW Files Version: 2.0 | November 2025 Use freely with attribution © K. R. Gunawardena Ph.D (Cantab), TU Delft, Delft, Netherlands. 1.0 Purpose and Scope This tool was initially developed as part of a project funded by the University of Bath (Version: 1.0, attribution should be made to Gunawardena et al. (2019, 2017)). The current version describes the revised methodology to compute hourly reference evapotranspiration (𝐸𝑇) for a short reference crop (grass) using the FAO-56 Penman-Monteith formulation (Eq. 53), and meteorological drivers extracted from EPW (EnergyPlus Weather) files. The method generates 𝐸𝑇 in mm·h⁻¹ and is suitable for urban microclimate and water-balance assessments where hourly resolution is required. 2.0 Theoretical background This section elaborates on the physical basis of the FAO-56 Penman–Monteith (PM) equation at an hourly time step. It complements Section 3 and keeps the same symbols and units used there (kPa, MJ·m⁻²·h⁻¹, mm·h⁻¹) (Allen et al. 1998). 2.1 Equations and conventions Hourly FAO-56 Penman-Monteith (Eq. 53): 𝑬𝑻𝟎=0.408 · 𝛥 · (𝑅𝑛 − 𝐺)+ 𝛾 · 37/(𝑇 + 273.15)· 𝑢₂ · (𝑒𝑠 − 𝑒𝑎) 𝛥 + 𝛾 · (1 + 0.34 · 𝑢₂) (mm · h⁻¹) Where: 𝑻 is air temperature (°C), 𝒖₂ is wind speed at 2 m (m·s⁻¹), 𝒆𝒔 and 𝑒𝑎 are saturation and actual vapour pressure (kPa), 𝜟 is the slope of the saturation vapour pressure curve (kPa·°C⁻¹), 𝜸 is the psychrometric constant (kPa·°C⁻¹), 𝑹𝒏 is net radiation (MJ·m⁻²·h⁻¹), and 𝑮 is soil heat flux (MJ·m⁻²·h⁻¹). 2.2 Reference surface and resistances FAO-56 defines the short reference crop as a hypothetical, well-watered grass canopy (height ≈ 0.12 m), with a fixed canopy (surface) resistance 𝑟𝑐 = 70 𝑠 · 𝑚⁻¹ and albedo 𝛼 = 0.23. These standardised parameters isolate the atmospheric demand component of 𝐸𝑇₀ and ensure results are comparable across sites and periods. In our workflow these values are configurable but default to the FAO recommendations. Page 2 of 6 2.3 Psychrometrics and air pressure Hourly saturation vapour pressure 𝑒𝑠(𝑇) and its slope 𝛥(𝑇) are evaluated from air temperature using the FAO-56 formulations. The psychrometric constant 𝛾 depends on the moist-air pressure 𝑃, which we estimate from the site elevation via the barometric formula; this captures the reduction of 𝛾 with altitude and is important in mountainous regions.  Saturation vapour pressure: 𝑒𝑠(𝑇)= 0.6108 · 𝑒𝑥𝑝[17.27 · 𝑇/(𝑇 + 237.3)] (kPa). Slope 𝛥(𝑇)= 4098 · 𝑒𝑠(𝑇)/(𝑇 + 237.3) (kPa·°C⁻¹).  Psychrometric constant: 𝛾 = 0.000665 · 𝑃, where 𝑃 (kPa) is atmospheric pressure estimated from elevation: 𝑃 = 101.3 · (293 − 0.0065 · 𝑧)/293.26. 2.4 Aerodynamics: wind adjustment to 2 m Because the PM equation uses wind speed at 2 m (𝑢₂), observed wind at the EPW anemometer height 𝑧ᵤ is converted via the FAO neutral log-profile relationship. When 𝑧ᵤ = 10 m (typical for EPW), 𝑢₂ ≈ 𝑢₁₀ × 4.87 / 𝑙𝑛(67.8 · 10 − 5.42). A minimum 𝑧ᵤ = 0.1 m is enforced to avoid singularities. 2.5 Radiation terms and hourly solar geometry Shortwave at the surface Rs is obtained from EPW global horizontal radiation (Wh·m⁻²) after unit conversion to MJ·m⁻²·h⁻¹. Net shortwave is 𝑅𝑛𝑠 = (1 − 𝛼) · 𝑅𝑠. Hourly extra-terrestrial radiation 𝑅𝑎 is computed from solar declination, inverse Earth–Sun distance, and the start/end hour angles; this yields a physically consistent clear-sky radiation 𝑅𝑠𝑜 = (0.75 + 2 × 10⁻⁵ · 𝑧) · 𝑅𝑎 that scales with altitude. Net longwave 𝑅𝑛𝑙 adapts the FAO daily emissivity-cloud term to the hourly step using σ/24 and a bounded cloud factor from 𝑅𝑠/𝑅𝑠𝑜; bounds avoid spuriously large magnitudes at night (𝑅𝑠𝑜 → 0).  Shortwave: 𝑅𝑠 = 𝐺𝐻𝐼 · 0.0036 (Wh·m⁻² → MJ·m⁻²·h⁻¹). Net shortwave: 𝑅𝑛𝑠 = (1 − 𝛼)· 𝑅𝑠, with surface albedo 𝛼.  Extraterrestrial radiation (hourly): 𝑅𝑎 computed from solar geometry (declination, inverse Earth–Sun distance, and hour angles at the start/end of each hour). Clear-sky radiation: 𝑅𝑠𝑜 = (0.75 + 2 × 10 · 𝑧)· 𝑅𝑎, with 𝑧 in m.  Net longwave: 𝑅𝑛𝑙 = 𝜎 · 𝑇 · 0.34 − 0.14 · √𝑒𝑎  · (1.35 · 𝑅𝑠/𝑅𝑠𝑜 − 0.35), using 𝜎 = 4.903 × 10 MJ·K⁻⁴·m⁻²·day⁻¹, adapted to an hourly constant ( 𝜎/24). Cloud term is bounded to avoid non-physical values at night. 2.6 Temporal treatment and soil heat flux EPW timestamps represent the end of each hour; for solar geometry calculations, the mid-hour (−30 min) is used to ensure correct sun angles. Soil heat lux (𝐺) is not computed from a full soil heat budget but approximated as a fixed fraction of net radiation (𝑅𝑛), following FAO-56 guidance: 𝐺 = 𝑓 · 𝑅𝑛, where 𝑓 = 0.10 during daylight and 𝑓 = 0.50 at night by default. These fractions are user-configurable and provide a simple way to capture diurnal phase lag without complex modeling. Page 3 of 6 3.0 Quality control and numerical safeguards These safeguards improve robustness when EPW fields are missing or in low-sun conditions.  Bound 𝑅𝑠/𝑅𝑠𝑜 within [0, 1.0] when 𝑅𝑠𝑜 > 0; set to 0, when 𝑅𝑠𝑜 = 0.  Bound cloud factor in 𝑅𝑛𝑙 to [0.05, 1.0], and clip 𝑅𝑛𝑙 ≥ 0.  Clip 𝑉𝑃𝐷 ≥ 0 and 𝐸𝑇≥ 0; replace non-finite values with 0.  Check presence of dewpoint; use 𝑅𝐻 fallback when not available.  Ensure anemometer height ≥ 0.1 m in wind conversion. 4.0 Data Inputs  EPW variables: Global Horizontal Radiation (Wh m⁻²), Dry-bulb Temperature (°C), Dewpoint (°C, preferred), Relative Humidity (%), Wind speed (m s⁻¹), Date-time fields.  Site parameters (Excel): latitude (deg), longitude (deg), time-zone UTC offset (h), elevation (m), surface albedo (–), anemometer height (m), soil heat flux fractions (day/night), dew-point usage flag. 5.0 Processing workflow 1. Parse EPW; create mid-hour timestamps (end-of-hour minus 30 minutes). 2. Compute 𝑅𝑠,𝑒𝑠(𝑇),𝑒𝑎 (from dewpoint if flagged, else from 𝑅𝐻 and 𝑇), 𝑉𝑃𝐷 = 𝑚𝑎𝑥(𝑒𝑠 − 𝑒𝑎,0). 3. Convert wind speed at measurement height to 2 m using FAO log profile. 4. Compute 𝛥(𝑇) and 𝛾(𝑧). 5. Derive solar geometry per hour; compute 𝑅𝑎 and 𝑅𝑠𝑜; then 𝑅𝑛𝑠,𝑅𝑛𝑙,𝑎𝑛𝑑 𝑅𝑛. 6. Estimate 𝐺 via day/night fractions; finally compute 𝐸𝑇 using Eq. 53. 7. Clip negative 𝐸𝑇 to 0 and write all diagnostics to Excel. 6.0 Execution workflows 6.1 MATLAB workflow  Function: PM56_ET0_hourly (epwFile, paramInputXLSX, "ResultFile", "ET0_hourly_results.xlsx").  Inputs: (a) EPW file (.epw or EPW-formatted .txt); (b) Excel “Parameters” sheet (data from row 4). If latitude/longitude/timezone/elevation are blank in Excel, the EPW LOCATION header is used.  Run: from MATLAB current folder or addpath to the function; specify the EPW and Excel files.  Output: an Excel workbook with three tabs: Results (𝐸𝑇 𝐻𝑜𝑢𝑟𝑙𝑦: hourly 𝐸𝑇 mm·h⁻¹, and diagnostics 𝑅𝑛,𝐺,𝑅𝑠,𝑅𝑠𝑜,𝑇,𝑢2,𝑒𝑠,𝑒𝑎,𝑉𝑃𝐷,𝛥,𝛾), Output Parameters (descriptions and units), CalcNotes (simulation notes). Page 4 of 6 6.2 Python workflow • Script: PM56_ET0_hourly.py (command-line tool). Dependencies: pandas, numpy, openpyxl. Install with: pip install pandas numpy openpyxl. • Inputs: same as MATLAB (EPW and Excel “Parameters” sheet starting at row 4). • Run: python PM56_ET0_hourly.py <EPW_FILE> <PARAMS_XLSX> --out ET0_hourly_results.xlsx • Output: Excel workbook with tabs: Results (𝐸𝑇 𝐻𝑜𝑢𝑟𝑙𝑦: hourly 𝐸𝑇 mm·h⁻¹, and diagnostics 𝑅𝑛,𝐺,𝑅𝑠,𝑅𝑠𝑜,𝑇,𝑢2,𝑒𝑠,𝑒𝑎,𝑉𝑃𝐷,𝛥,𝛾 ), Output Parameters (descriptions and units), CalcNotes (simulation notes). 7.0 Assumptions and limitations Assumes reference short crop (grass) with 𝑟𝑐 = 70 s·m⁻¹ and albedo 0.23. Hourly net longwave formula adapts FAO daily cloud factor; results are sensitive to radiation quality and to the dewpoint vs 𝑅𝐻 choice. EPW wind is assumed at the given anemometer height; if unknown, 10 m is typical. 8.0 MATLAB vs. Python verification test In this sample Amsterdam IWEC.epw file the nearest mid-summer date with a complete day is 1996-06-21. For this date, the MATLAB column is computed with the MATLAB algorithm reproduced in Python (parity check). Date & Time (end of hour) 𝑬 𝑻 𝟎 (Python) [mm h⁻¹] 𝑬 𝑻 𝟎 (MATLAB) [mm h⁻¹] Abs Diff [mm h⁻¹] Rel Diff [%] 1996-06-21 00:00:00 0.033204 0.033204 0.000000 0.000000 1996-06-21 01:00:00 0.035171 0.035171 0.000000 0.000000 1996-06-21 02:00:00 0.040081 0.040081 0.000000 0.000000 1996-06-21 03:00:00 0.028787 0.028787 0.000000 0.000000 1996-06-21 04:00:00 0.049238 0.049238 0.000000 0.000000 1996-06-21 05:00:00 0.046134 0.046134 0.000000 0.000000 1996-06-21 06:00:00 0.045868 0.045868 0.000000 0.000000 1996-06-21 07:00:00 0.071889 0.071889 0.000000 0.000000 1996-06-21 08:00:00 0.112356 0.112356 0.000000 0.000000 1996-06-21 09:00:00 0.167223 0.167223 0.000000 0.000000 1996-06-21 10:00:00 0.205673 0.205673 0.000000 0.000000 1996-06-21 11:00:00 0.157285 0.157285 0.000000 0.000000 1996-06-21 12:00:00 0.191063 0.191063 0.000000 0.000000 1996-06-21 13:00:00 0.220719 0.220719 0.000000 0.000000 1996-06-21 14:00:00 0.192751 0.192751 0.000000 0.000000 1996-06-21 15:00:00 0.150901 0.150901 0.000000 0.000000 1996-06-21 16:00:00 0.187324 0.187324 0.000000 0.000000 1996-06-21 17:00:00 0.184559 0.184559 0.000000 0.000000 1996-06-21 18:00:00 0.110979 0.110979 0.000000 0.000000 1996 - 06 - 21 19:00:00 0.116553 0.116553 0.000000 0.000000 1996-06-21 20:00:00 0.085613 0.085613 0.000000 0.000000 1996-06-21 21:00:00 0.035523 0.035523 0.000000 0.000000 1996-06-21 22:00:00 0.027081 0.027081 0.000000 0.000000 1996-06-21 23:00:00 0.024874 0.024874 0.000000 0.000000 Page 5 of 6 9.0 Potential uses of this model • Urban microclimate assessment: quantify atmospheric water demand for thermal comfort and heat-mitigation studies. • Irrigation scheduling and water-balance modelling for urban greenspaces, parks, and green roofs (hourly granularity). • Design and evaluation of nature-based solutions (street trees, bioswales, rain gardens) under different climates. • Scenario analysis for climate resilience: sensitivity of 𝐸𝑇₀ to heatwaves, humidity anomalies, and wind regimes. • Hydrological modelling linkages: coupling 𝐸𝑇₀ with soil-plant-atmosphere or urban drainage models at sub-daily steps. • Data QA/QC: use drivers and diagnostics (𝑅𝑛,𝐺,𝑉𝑃𝐷,𝑢₂) to screen EPW records and identify inconsistencies. 10.0 Open sharing license with attribution Copyright © 2025 Kanchane Gunawardena Permission is hereby granted, free of charge, to any person obtaining a copy of this work and associated documentation files (the “Work”), to use, reproduce, modify, publish, distribute, and create derivative works of the Work, subject to the following conditions: 1. Attribution You must give appropriate credit, provide a link to the license, and indicate if changes were made. Attribution must be clearly visible in any public use, publication, or redistribution of the Work. Suggested citation format: “Hourly FAO-56 Penman–Monteith Reference Evapotranspiration (𝐸𝑇) from EPW Files, by K. Gunawardena, licensed under the Open Sharing License with Attribution.” 2. No Warranty The Work is provided “as is”, without warranty of any kind, express or implied, including but not limited to the warranties of merchantability, fitness for a particular purpose, and noninfringement. 3. Redistribution Copies or substantial portions of the Work must retain this license and attribution notice. 4. Modifications You may modify the Work, but must clearly indicate that changes were made and retain attribution to the original author. 11.0 References Allen, R.G., Walter, I.A., Elliott, R.L. et al. (2005) The ASCE Standardized Reference Evapotranspiration Equation. ASCE-EWRI Task Committee Report, ASCE. Allen, R.G., Pruitt, W.O., Wright, J.L. et al. (2006) A recommendation on standardized surface resistance for hourly calculation of reference ET0 by the FAO-56 Penman–Monteith method. Agricultural Water Management, 81(1–2), 1–22. Page 6 of 6 Allen, R. G., Pereira, L. S., Raes, D., & Smith, M. (1998). Crop evapotranspiration-Guidelines for computing crop water requirements-FAO Irrigation and drainage paper 56. FAO, Rome, 300(9), D05109. Gunawardena, K., Kershaw, T., & Steemers, K. (2019). Simulation pathway for estimating heat island influence on urban/suburban building space-conditioning loads and response to facade material changes. Building and Environment, 150(January), 195–205. Gunawardena, K. R., Mccullen, N., & Kershaw, T. (2017). Heat island influence on space-conditioning loads of urban and suburban office buildings. In Cities and Climate Conference 2017, Potsdam: Potsdam Institute for Climate Impact Research, pp. 1–13.