The study aimed to evaluate the accuracy of empirical equations (Hargreaves-Samani; HS, Irmak; IR and Dalton; DT) and multivariate linear regression models (MLR1–6) for estimating reference evapotranspiration (ETRef) in different climates of Iran based on the Köppen method including arid desert (Bw), semiarid (Bs), humid with mild winters (C), and humid with severe winters (D). For this purpose, climatic data of 33 meteorological stations during 30 statistical years 1990–2019 were used with a monthly time step. Based on various meteorological data (minimum and maximum temperature, relative humidity, wind speed, solar radiation, extraterrestrial radiation, and vapor pressure deficit), in addition to 6 multivariate linear regression models and three empirical equations were used as MLR1, MLR2, and HS (temperature-based), MLR3 and IR (radiation-based), MLR4, MLR5 and DT (mass transfer-based), and MLR6 (combination-based) were also used to estimate the reference evapotranspiration. The results of these models were compared using the root mean square error (RMSE), mean absolute error (MAE), scatter index (SI), determination coefficient (R2), and Nash-Sutcliffe efficiency (NSE) statistical criteria with the evapotranspiration results of the FAO56 Penman-Monteith reference as target data. All MLR models gave better results than empirical equations. The results showed that the simplest regression model (MLR1) based on the minimum and maximum temperature data was more accurate than the empirical equations. The lowest and highest accuracy related to the MLR6 model and HS empirical equation with RMSE was 10.8–15.1 mm month−1 and 22–28.3 mm month−1, respectively. Also, among all the evaluated equations, radiation-based models such as IR in Bw and Bs climates with MAE = 8.01–11.2 mm month−1 had higher accuracy than C and D climates with MAE = 13.44–14.48 mm month−1. In general, the results showed that the ability of regression models was excellent in all climates from Bw to D based on SI < 0.2.