Being one of the most spectacular events in nature, lightning is basically a transient high current electric discharge in the atmosphere, which extends up to kilometres. Cloud to ground discharge is the most hazardous one as far as ground based structures are considered. Among the different phases of a lightning flash, return stroke is considered to be the most energetic phase and is basically responsible for most of the damages. Hence, much emphasis has been given to return stroke modeling. A more realistic modeling of return stroke is very essential to accurately study the interaction of return stroke with the structures on ground. As return stroke is dominated by electromagnetic phenomenon, an electromagnetic model will be the most suitable one. It does not call for any assumption on the mode of wave propagation, as well as, electromagnetic coupling between the different channel portions. There are mainly two approaches adopted for electromagnetic models i.e. frequency domain and time domain approach. Time domain approach is more reliable as it can handle, in principle, the nonlinear processes in the lightning channel. It is also free of numerical frequency domain to time domain inversion problem, which are found to be quite severe. However, most of the previous works on time domain electromagnetic models suffered from following two serious limitations - (i) the initial charge on the channel, which forms the true excitation for the problem, is not considered and (ii) instead of the non-linearly rising conductivity of the channel, a constant resistance is employed. For a realistic simulation of the interaction between the channel and any intercepting system, a time domain model with the above two major aspects being fully represented is very essential. In an earlier work, all these aspects have been fully considered but a domain based numerical modelling was employed. Consequently, it was difficult to consider the down conductor and further the number of unknowns was considerably large. In view of this, the present work is taken up and its scope is defined as to develop a boundary based numerical time domain electromagnetic model in which the initial charge on the channel and the non-linearly evolving channel conductance are fully considered. For the electrical engineering applications, electromagnetic aspects of the lightning phenomena is more important than the other associated physical processes and hence, importance is given only to the electromagnetic aspects. In other words, the light emission, thunder, chemical reactions at the channel etc. are not considered. Also, for most of the electrical engineering applications, the critical portion of current would be the region up to and around the peak and hence, modeling for this regime will be given prime importance. Owing to the complexity of the problem, some simplifying assumptions would be very essential. The literature indicates that these assumptions do not affect the adequate representation of the phenomena. Lightning channel is considered to be vertically straight without any branches. Earth is considered to be perfectly conducting. Explicit reference to dynamically varying channel radius, temperature and the air density is not made. However, it is assumed that the arc equation employed to describe the temporal changes in conductivity would adequately take care of these parameters. Lightning channel is represented by a highly conducting small core, which is surrounded by a weakly conducting corona sheath. The initial charge on the channel is deduced by solving for electrostatic field, with leader portion set to possess an axial gradient of 60 V/cm and the streamer portion to 5 – 10 kV/cm. The radius of the corona sheath is set iteratively by enforcing a gradient of 24 kV/cm up to its radial boundary. As analytical solution for the problem is impractical, suitable numerical solution is sought. Since the spatial extension of this time marching problem is virtually unbound and that the significant conduction is rather solely confined to an extremely small cross section of the channel core, a boundary-based method is selected. Amongst the numerical methods, the present work employs the moment method for the solution of the fields associated with the return strokes. A numerical solution of the Electric Field Integral Equation (EFIE) for thin structures has been developed in the literature. The same approach has been employed in the present work, however, with suitable modifications to suit the lightning problem. The code was written in MATLAB and integrations involved in the EFIE were solved using MATLAB symbolic computation. Before introducing the channel dynamic conductance and the initial charge on the channel, the code developed is validated by comparing the results for a center fed dipole antenna with that given in the literature. Also, NEC (Numeric Electromagnetic Code) simulations for various cases of monopole and dipole antenna were performed. The results from the code developed are shown to have good matching with that obtained from NEC based time domain results. In an earlier work, the dynamic conductance of the return stroke channel core, which is a high current electric arc, was represented by a first order arc equation. The same approach is employed in the present work also. Similarly, the transition from streamer to leader was modeled by Braginskii’s spark law and the same has been considered in the present work. A value of 10-5 S/m was used for minimum value of streamer conductance. For numerical stability, upper (Gmax = 3 S/m) and lower bounds (Gmin = 0.0083 S/m) for the channel conductance are forced. Preliminary simulations were run with and without dynamic channel conductance. The initial charge distribution along the channel formed the excitation. Results clearly show that without the dynamically varying channel conductance, no streamer to leader transition and hence, no return stroke evolution can occur. In other words, the non-linearly evolving channel conductance is mainly responsible for the evolution of the return stroke. In order to consider the charge neutralization by the return stroke, the charge deposited by it is diffused into the corona sheath. A fixed value of the corona sheath conductance is employed and the diffusion process is modeled by an equation derived from the continuity equation. To study the effect of corona sheath, simulations were run with and without corona. From the simulation results it was observed that the corona sheath causes increase in peak value of the stroke current, as well as, time to front and a decrease in the velocity of propagation. For the validation of the model, the basic characteristics of the return stroke current like the current wave shape, temporal variation of stroke current at different heights, velocity of propagation and the vertical electric fields at various radial distances were compared with available field/experimental data. A good agreement was seen and based on this, it is concluded that the present work has successfully developed a boundary based time domain numerical model for the lightning return stroke. Natural lightning being a stochastic process, the values of the parameters associated with it would differ in every event. On other hand, any deterministic model like the one developed in the present work predicts a fixed pattern of the simulated quantities. Therefore, it was felt that some of the model parameters must be permitted to vary so that a range of results could be obtained rather than a single pattern of results. Incidentally, the model parameters like arc time constant, settling value of arc conductivity/gradient, bounds for channel conductivity, streamer gradient, radius of the core etc. are not precisely known for the natural lightning environment. Further, some of them are known to vary within an event. Considering these and that simplicity is very important in already complex model, the above-mentioned parameters are taken as tunable parameters (of course to be varied within the prescribed range) for deducing the return stroke currents with some desired characteristics. A study on the influence of these parameters is made and suggestions are provided. Simulations for the nominal range of stroke currents are made and results are presented. These simulations clearly show the role of cloud potential, which in turn dictates the length of final bridging streamer, on the return stroke currents. The spatio-temporal variation of the current, charge deposited by the return stroke and the channel conductivity are presented which, reveal the dynamic processes leading to the evolution of return stroke current. Subsequently, simulations for two cases of stroke to elevated strike object are attempted. The upward leader was modeled quite similar to the descending one. Many interesting findings are made. In summary, the present work has successfully developed a boundary-based time domain numerical electromagnetic model for the lightning return stroke, wherein, the initial charge deposited on the channel and the non-linearly rising channel conductance have been appropriately considered. Simulation using the model clearly depicts the dynamic evolution of the return stroke. The characteristics of the simulated return strokes are in good agreement with the field data. Some of the parameters of the model are suggested as tunable parameters, which permit simulation of strokes with different characteristics.