What would be the effect of repeating samples in time domain on frequency domain?
It's quite benefical to view the mathematical expression of the process as well.

In discrete-time interpolation, a mere $N$ times repetition of a single sample (pixel) into the new enlarged block is an interpolation with zero-order hold kernel. This is mostly applied for simple fast but crude image zooming.

In 1D, the kernel is $$h_0[n] = \begin{cases}{ 1 ~~~,~~~ 0 \leq n \leq N-1 \\ 0 ~~~,~~~ \text{ otherwise } } \end{cases} \tag{1}$$

The interpolated signal can be expressed as :

$$y[n] = \sum_{k} x[k] h_0[n-kN] \tag{2}$$ this can be expressed as a convolution like 

$$y[n] = \left(\sum_{k} x[k] \delta[n- Nk] \right) \star h_0[n] \tag{3}$$

where the signal in parenthesis is the expanded version of $x[n]$:

$$x[n] \rightarrow \boxed{ \uparrow N} \rightarrow x_e[n] = \sum_{k} x[k] \delta[n- Nk] \tag{4}$$

Applying DT Fourier transform on Eq(3) yields:

$$ Y(e^{j\omega}) = X_e(e^{j\omega}) \cdot H_0(e^{j\omega})  \tag{5} $$

where $X_e(e^{j\omega}) = X(e^{j\omega N}) $ due to the expansion operation in Eq(4) and $H_0(e^{j\omega})$ is the DTFT of the rectangular pulse $h_0[n]$ which is:

$$H_0(w) = e^{-j\omega (N-1)/2} \frac{ \sin( \omega N /2 ) } { \sin( \omega / 2 ) } $$

then the result becomes:

$$ Y(e^{j\omega}) = X(e^{j\omega N}) \cdot e^{-j\omega (N-1)/2} \frac{ \sin( \omega N /2 ) } {\sin( \omega / 2 ) }  \tag{6} $$

Whose magnitude is :
$$ \boxed{ |Y(e^{j\omega})| = |X(e^{j\omega N})| \cdot \left| \frac{ \sin( \omega N /2 ) } {\sin( \omega / 2 ) } \right| }  \tag{7} $$

Note that the $X(e^{j \omega N})$ is composed of repeated images of $X(e^{j \omega})$. And the effect of the interpolation kernel $h_0[n]$ is to remove the images as much as possible with that crude frequency response $H_0(e^{j \omega}) $.

The following Octave/Matlab code demonstrates the operation:

clear all; close all; clc;

N = 8;    % interpolation
L = 32;   % original signal length
x = 5*( sin(0.3*pi*([0:L-1]-L/2))./(pi*([0:L-1]-L/2))).^2;
x(isnan(x))=5*0.3*0.3;

h0 = ones(1,N);

xe = zeros(N*L,1);  % expanded signal
xe(1:N:N*L) =  x;

yi = conv(xe, h0);  % the interpolated signal: y[n] = xe[n] * h0[n]

figure,subplot(2,1,1);
stem(x);
title('original signal');
subplot(2,1,2);
stem(yi);
title('zero-order hold interpolated output');


figure,subplot(2,1,1);
plot(linspace(0,2,1024),abs(fft(x,1024)),'b');
title('Spectrum magnitude of x[n]')

subplot(2,1,2);
plot(linspace(0,2,1024),abs(fft(xe,1024)),'b');
hold on
plot(linspace(0,2,1024),abs(fft(h0,1024)),'r');
plot(linspace(0,2,1024),abs(fft(yi,1024)),'g','linewidth',2);
legend('|X_e(w)|','|H_0(w)|','|Y_i(w)|');


with the following outputs:



and the spectrum:



as you can see, the nulls of the periodic sinc (the Dirichlet kernel) are at the centers of the repeated spectra (the images) of $X(w)$. However this is not enough to suppress them all unlike an ideal lowpass filter (or one of its adequate approximations) which would completely eliminate the image spectra and yield only the one inside the main lobe; hence yielding a perfect interpolation compared to this crude zero order hold kernel giving.