The distribution of slip during an earthquake and how it propagates among faults in the subduction system play a major role in seismic and tsunami hazards, yet they are poorly understood because offshore observations are often lacking. Here we derive the slip distribution and rupture evolution during the 2016 M-w 7.9 Kaikoura (New Zealand) earthquake that reconcile the surface rupture, space geodetic measurements, seismological and tsunami waveform records. We use twelve fault segments, with eleven in the crust and one on the megathrust interface, to model the geodetic data and match the major features of the complex surface ruptures. Our modeling result indicates that a large portion of the moment is distributed on the subduction interface, making a significant contribution to the far field surface deformation and teleseismic body waves. The inclusion of local strong motion and teleseismic waveform data in the joint inversion reveals a unilateral rupture towards northeast with a relatively low averaged rupture speed of similar to 1.5 km/s. The first 30 s of the rupture took place on the crustal faults with oblique slip motion and jumped between fault segments that have large differences in strike and dip. The peak moment release occurred at similar to 65 s, corresponding to simultaneous rupture of both plate interface and the overlying splay faults with rake angle changes progressively from thrust to strike-slip. The slip on the Papatea fault produced more than 2 m of offshore uplift, making a major contribution to the tsunami at the Kaikoura station, while the northeastern end of the rupture can explain the main features at the Wellington station. Our inversions and simulations illuminate complex up-dip rupture behavior that should be taken into consideration in both seismic and tsunami hazard assessment. The extreme complex rupture behavior also brings new challenges to the earthquake dynamic simulations and understanding the physics of earthquakes.