## 1 Introduction  
Irreversible Electroporation (IRE) is a local ablation modality that is currently investigated for the treatment of unresectable tumors, including liver, lung, pancreatic, and urological cancers [1; 2; 3; 4; 5]. Contrary to conventional thermal ablation techniques, in which ablation is achieved by exceeding the thermal-damage temperature threshold, the IRE procedure is based on the application of short, high-voltage pulses across needle or plate electrode pairs implanted in and surrounding the target volume to cause nanopores in the cell membranes [6]. The pore formation induces disturbances in the cell's homeostasis, with as a consequence cell death through accidental and regulated cell death mechanisms, and pore formation also improves solute mobility which translates in increased electrical conductivity in the bulk tissue [7; 8]. IRE ablation depends on several parameters, such as the electrical pulse parameters (voltage, shape, duration, number, and frequency), temperature, and tissue properties of the target volume [9; 10; 11; 12]. Therefore, it is important to optimize the applied IRE protocols to mainly obtain ablation through electroporation only or the combined effect of electroporation and mild hyperthermia, without achieving thermal damage [12]. In addition to pre-clinical in vitro and in vivo experiments, computational modeling provides a powerful tool for investigating and predicting the outcome of IRE-protocols. Using computational models, a wide variety of IRE protocol variations can be explored. Specifically, the influence of individual parameters can easily be investigated by varying a single parameter, while maintaining other parameters constant. This way, computational modelling can assist in the optimization of IRE protocols. Currently, several commercial and non-commercial software packages are available to calculate distributions of physical key quantities that occur during an electroporation treatment, such as electric field and temperature, and the change in behavior of tissue properties as a function of the former quantities [12; 13]. Examples of packages capable of implementing field and temperature dependent tissues properties include COMSOL Multiphysics (COMSOL AB, Stockholm, Sweden), QuickField (Terra Analysis, Svendborg, Denmark), Eview ([[https://eview.upf.edu/](https://eview.upf.edu/)]([https://eview.upf.edu/](https://eview.upf.edu/)); accessed on 4 May 2021) [14], IRENA ([[https://team.inria.fr/monc/software/](https://team.inria.fr/monc/software/)]([https://team.inria.fr/monc/software/](https://team.inria.fr/monc/software/)); accessed on 4 May 2021) [15], and OpenEP ([[https://github.com/LSC-UBA/OpenEP/](https://github.com/LSC-UBA/OpenEP/)]([https://github.com/LSC-UBA/OpenEP/](https://github.com/LSC-UBA/OpenEP/)); accessed on 28 February 2021) [13]. Depending on the software package, the numerical simulations are generally performed using the Finite Element Method and/or the Finite Difference Method (FDM). Some software packages allow for local adaptive mesh refinement in the vicinity of small or round structures to improve the modelling accuracy of needle electrodes or cylindrically shaped electrodes [16; 17]. However, a uniform rectangular grid is usually used in FDM-based software packages [13; 15; 18; 19]. The small diameter of the needle electrode applied in the clinic may result in reduced accuracy of the calculated electric potential in the vicinity of the needle electrode model in the case of an inappropriate choice of mesh properties in FDM-based software packages. Particularly, this is the case when large grid voxels (from here on the term \"voxels\" will be used instead of \"mesh elements\" as a reference for the elements of a uniform rectangular grid) are obtained from Computed Tomography (CT) scans without increasing the grid resolution.